Seeking advice on saltwater intrusion modelling with unstructured meshes

[Micaela Raviola]

Dear all,

I am currently working on reactive transport modelling in coastal aquifers.

To simulate the dynamics of saline wedge intrusion, I need to consider density-driven flow. I am currently using the AUXILIARY_SALINITY keyword to calculate density-dependent flow. However, as indicated in the documentation, AUXILIARY_SALINITY is incompatible with the HYDROSTATIC flow mode.
To overcome this, I am using a Dirichlet boundary condition, calculating the pressure for each boundary cell externally and importing it into the simulation as a gridded dataset (similar to the multiphase_regional_flow example approach). The main challenge lies in replicating this approach with unstructured meshes.
As gridded datasets cannot be used for unstructured meshes and cell-indexed datasets are not supported for boundary conditions, I am looking for the best way to map variable pressures as boundary conditions for an unstructured mesh domain.

How can I apply spatially variable boundary conditions (calculated externally) to an unstructured mesh while continuing to use AUXILIARY_SALINITY?

Alternatively, is there a different workflow or flow mode for simulating saltwater intrusion in unstructured grids that would avoid the conflict between AUXILIARY_SALINITY and HYDROSTATIC, or between cell-indexed datasets and boundary conditions?

Any advice, examples or alternative solutions would be greatly appreciated.

Thank you in advance. 

Micaela Raviola

[Hammond, Glenn E]

Micaela,

Have you tried to experiment with gridded datasets on unstructured grids? If so, did you encounter any error messages during your attempts? The key advantage of these datasets lies in their uniform, structured format, which allows the code to perform linear interpolation of values to the face centers—even when working within unstructured grid environments.

Additionally, I suggest running a 1D vertical column simulation with salinity, uniform grid spacing, and a prescribed boundary pressure until it reaches steady state (long time). You can then use these results in a gridded dataset, much like the multiphase regional double approach. This method provides a more accurate hydrostatic profile, since fluid density varies with elevation and influences the pressure gradient.

Please let me know if this approach is helpful.

Glenn

From: pflotra...@googlegroups.com <pflotra...@googlegroups.com> on behalf of Micaela Raviola <micaela...@gmail.com>
Date: Wednesday, March 11, 2026 at 7:47 AM
To: pflotran-users <pflotra...@googlegroups.com>
Subject: [pflotran-users: 8705] Seeking advice on saltwater intrusion modelling with unstructured meshes

You don't often get email from micaela...@gmail.com. Learn why this is important

Check twice before you click! This email originated from outside PNNL.

--
You received this message because you are subscribed to the Google Groups "pflotran-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pflotran-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/pflotran-users/9d45a0cd-31fd-4fc1-a909-0e4df574ddaen%40googlegroups.com.

[Micaela Raviola]

Glenn,
Thank you very much for your suggestions. I have just tried using gridded datasets on unstructured meshes and it works perfectly.
Just to make sure I understand correctly, are you suggesting running a 1D simulation over a long period of time while already considering salinity? For example, imposing a Dirichlet boundary condition on the bottom (considering the density of seawater), and then using the steady state result as input for the aquifer simulation?

 Thank you again!

Micaela

[Hammond, Glenn E]

That’s correct. Begin by creating a one-dimensional column, specifying the salinity using your previous approach. Next, prescribe a pressure at either the top or bottom of the column to reflect the desired field conditions.

Are you looking for the boundary condition to vary over time?

Glenn

To view this discussion visit https://groups.google.com/d/msgid/pflotran-users/76216aef-d7e9-49cd-98d7-5a854bf4dbdan%40googlegroups.com.

[Micaela Raviola]

Yes, I’d definitely like to try simulating time-varying boundary conditions soon. At the moment, I’m trying to achieve a result that’s consistent with the physics of the problem. I’m having trouble optimising the parameters related to maximum timestep size  and the numerical methods since the flow is struggling to converge. Do you have any suggestions on how to tune these parameters most effectively?

Once I’ve built a working infrastructure, I’d like to implement the time-varying conditions. Do you have any suggestions for this? 

Thank you so much for your previous insights. 

Micaela

[Hammond, Glenn E]

Micaela,

Could you please send the screen output from your simulation? You can capture this by running:

/path/to/pflotran -input_prefix <prefix> | tee pflotran.stdout

If the output is too extensive, feel free to share just the timesteps that are having trouble converging.

Glenn

To view this discussion visit https://groups.google.com/d/msgid/pflotran-users/7246456a-8dc4-4f3b-a05c-000cd9f54909n%40googlegroups.com.

[Micaela Raviola]

Glenn,

 here are the kinds of convergence issues that periodically pop up during the simulation. I’m only attaching two screenshots here because the output is quite long.

 Initially, I thought it was a problem related to the initial time step of the simulation. However, I realised that it keeps happening (periodically) long after the simulation has started (even though the simulation has come to an end).

Thanks again for your time!

Micaela 

== GENERAL MULTIPHASE FLOW =====================================================
  0 2r: 5.90E-05 2x: 0.00E+00 2u: 0.00E+00 ir: 7.84E-06 iu: 0.00E+00 rsn:   0
  1 2r: 1.67E-05 2x: 8.41E+06 2u: 5.26E+01 ir: 1.28E-05 iu: 5.17E+00 rsn:   0
  2 2r: 1.20E-04 2x: 8.41E+06 2u: 2.09E+00 ir: 9.18E-05 iu: 2.09E+00 rsn:   0
  3 2r: 3.66E-05 2x: 8.41E+06 2u: 4.55E+00 ir: 2.81E-05 iu: 4.55E+00 rsn:   0
  4 2r: 1.20E-04 2x: 8.41E+06 2u: 4.55E+00 ir: 9.18E-05 iu: 4.55E+00 rsn:   0
  5 2r: 3.66E-05 2x: 8.41E+06 2u: 4.55E+00 ir: 2.81E-05 iu: 4.55E+00 rsn:   0
  6 2r: 1.20E-04 2x: 8.41E+06 2u: 4.55E+00 ir: 9.18E-05 iu: 4.55E+00 rsn:   0
  7 2r: 3.66E-05 2x: 8.41E+06 2u: 4.55E+00 ir: 2.81E-05 iu: 4.55E+00 rsn:   0
 WARNING: Potential oscillatory convergence
  8 2r: 1.20E-04 2x: 8.41E+06 2u: 4.55E+00 ir: 9.18E-05 iu: 4.55E+00 rsn: -88
 WARNING: Potential oscillatory convergence
 -> Cut time step: snes= -88 icut= 1 [4] t= 3.74331E-04 dt= 3.79187E-06 [y]
Newton solver reason: Unknown(-88).
  0 2r: 5.68E-05 2x: 0.00E+00 2u: 0.00E+00 ir: 5.73E-06 iu: 0.00E+00 rsn:   0
  1 2r: 1.55E-06 2x: 8.41E+06 2u: 1.91E+01 ir: 1.19E-06 iu: 2.12E+00 rsn:   0
  2 2r: 3.39E-08 2x: 8.41E+06 2u: 4.29E-02 ir: 2.61E-08 iu: 2.86E-02 rsn: 999

 Step     14 Time=  3.78123E-04 Dt=  3.79187E-06 [y] conv_reason: 999
  newton =  10 [      65] linear =   262 [      2187] cuts =  1 [   4]
  --> SNES Linear/Non-Linear Iterations = 76 / 2
  --> SNES Residual: 0.339481E-07 0.113160E-10 0.260590E-07
  --> max chng: dpl=   3.7598E+00 dpg=   2.1190E+00 dpa=   2.1208E+00
                dxa=   8.6763E-09  dt=   1.3345E-04 dsg=   5.0649E-05

== GENERAL MULTIPHASE FLOW =====================================================
  0 2r: 6.45E-07 2x: 0.00E+00 2u: 0.00E+00 ir: 3.45E-07 iu: 0.00E+00 rsn:   0
  1 2r: 2.93E-09 2x: 8.41E+06 2u: 2.32E+00 ir: 1.50E-09 iu: 1.42E+00 rsn:   0
  2 2r: 6.84E-07 2x: 8.41E+06 2u: 3.94E+00 ir: 4.25E-07 iu: 2.76E+00 rsn:   0
  3 2r: 1.55E-08 2x: 8.41E+06 2u: 3.14E+00 ir: 8.15E-09 iu: 2.75E+00 rsn:   0
  4 2r: 1.14E-08 2x: 8.41E+06 2u: 2.41E+00 ir: 6.35E-09 iu: 1.84E+00 rsn:   0
  5 2r: 9.01E-08 2x: 8.41E+06 2u: 3.79E+00 ir: 5.64E-08 iu: 2.12E+00 rsn:   0
  6 2r: 5.23E-09 2x: 8.41E+06 2u: 3.37E+00 ir: 2.60E-09 iu: 1.40E+00 rsn:   0
  7 2r: 1.18E-08 2x: 8.41E+06 2u: 2.51E+00 ir: 6.64E-09 iu: 1.58E+00 rsn:   0
  8 2r: 4.83E-09 2x: 8.41E+06 2u: 1.96E+00 ir: 3.05E-09 iu: 1.57E+00 rsn: -88
 -> Cut time step: snes= -88 icut= 1 [112] t= 1.03370E+02 dt= 2.82919E-02 [y]
Newton solver reason: Unknown(-88).
  0 2r: 6.45E-07 2x: 0.00E+00 2u: 0.00E+00 ir: 3.45E-07 iu: 0.00E+00 rsn:   0
  1 2r: 1.35E-09 2x: 8.41E+06 2u: 1.40E+00 ir: 5.89E-10 iu: 8.18E-01 rsn:   0
  2 2r: 3.54E-07 2x: 8.41E+06 2u: 1.55E+00 ir: 2.21E-07 iu: 8.16E-01 rsn: 999

 Step   4905 Time=  1.03398E+02 Dt=  2.82919E-02 [y] conv_reason: 999
  newton =  10 [   12591] linear =   437 [    607735] cuts =  1 [ 112]
  --> SNES Linear/Non-Linear Iterations = 102 / 2
  --> SNES Residual: 0.353885E-06 0.117962E-09 0.220516E-06
  --> max chng: dpl=   1.1986E-01 dpg=   3.1780E-01 dpa=   3.1780E-01
                dxa=   1.8019E-08  dt=   4.6419E-04 dsg=   6.3512E-06

[Hammond, Glenn E]

Micaela,

Because you are operating in GENERAL mode, you can directly simulate salinity tightly coupled with flow. There is no need to use AUXILIARY_SALINITY, as it relies on loose coupling and is less effective than the tightly coupled method. However, you will still need to define a vertical hydrostatic boundary condition manually, since the code cannot generate this internally. While GENERAL mode is more complex than RICHARDS—the flow mode I originally assumed you were using—it remains a workable—perhaps better—option.

Glenn

To view this discussion visit https://groups.google.com/d/msgid/pflotran-users/b01c2bed-9a28-43a8-8e4e-da996d922e43n%40googlegroups.com.

[Micaela Raviola]

Hi Glenn,

I’ve managed to create a reasonable model (using GENERAL) of my simplified aquifer by leveraging gridded datasets and 1D column simulations to impose boundary conditions. However, I’m now trying to run the same simulation on a more complex, heterogeneous domain, and I’m encountering a significant issue that I can’t currently resolve. The 1D simulations that I use to stabilise the water table and then import as boundary conditions for my 3D domain are performed on material with homogeneous porosity and permeability. As long as the final 3D domain has the same homogeneous permeability as that used in the 1D simulations, everything works fine. However, when I try to impose the same boundary conditions on a heterogeneous domain, as shown in the attached image, I encounter convergence issues because the pressures I am imposing are inconsistent with the permeability distribution in the final domain.

Do you have any suggestions on how to work around this? For example, would it be effective to extract 1D columns from the final domain with the corresponding heterogeneous property distributions, or would you recommend an alternative approach?

Micaela