PRT fails to account for DRAPE option in MF6 model using Newton Under Relax Option

[javgs-bd]

We have been trying to implement PRT in a transient problem where the phreatic surface starts within a lower, relatively impermeable hydrogeological unit, and will raise in time to reach more permeable units. This is caused by a local recharge from the surface that starts at some time within the model simulation. We use the newton-raphson options to deal with the conversion of layers from dry to partially saturated. We have found that PRT crashes when releasing particles in the initially dry layer, regardles of the DRAPE option, when it should i) terminate the particles instead with ISTATUS 9 according to the docs (particle terminated immediately upon release into a dry cell) or ii) move the particles to the uppermost layer with partial saturation.

Describe the bug When using newton-raphson in the flow model, and particles are released in an initially dry layer, PRT crashes regardless of the DRAPE setting.

To reproduce In the attached notebook, ex-prt-tr.zip

1.- Go to the PRP section an uncomment this line, which makes particles to be released in layer 0. The uppermost and initially dry layer. # offsets = [(-1,0,0),(-1,-1,0), (-1,1,0), (-1,0,-1), (-1,0,1)] # particles are released in layer 0

The PRT simulation should crash with the message in screenshot 1, even when DRAPE is set to True.

2. In the GWF SIM and GWF Model section comment the following line. To deactivate NR option. newtonoptions = 'NEWTON UNDER_RELAXATION',

The PRT simulation would run to completion, particle trajectories change, but most importantly the flow solution is not the expected one as the upper layer remains dry. DRAPE seems to have been taken effect here.

3. An additional check can be made of activating back the NR option and setting DRAPE to False in the PRP section. The PRT simulation fails again with the same message (screenshot 2),

Expected behavior Particles should be either terminated or moved to the uppermost active layer according to DRAPE

Screenshots

Environment

• Operating system Windows 11
• MODFLOW 6 version 6.5.0 (screenshot 3)
  • Python 3.10.11 (tags/v3.10.11:7d4cc5a, Apr 5 2023, 00:38:17) [MSC v.1929 64 bit (AMD64)]
  • flopy version: 3.8.0

[wpbonelli]

Thanks for reporting this @javgs-bd. There are at least three separate issues here.

1) The div-0 FPE you are seeing has been fixed by a previous patch. The develop branch / nightly build is not affected by it.

2) There was a bug in the release time selection logic which produced additional/erroneous release times. This is fixed by #2017 which I will merge shortly.

However there is still a crash when newton is enabled and the particles are released from layer 0 (dry). With a layer 0 release, newton off, and drape on, the pathlines look reasonable; with drape off, the particles terminate immediately as expected. So this seems newton-specific. Currently, expected behavior with newton is to pass particles in dry cells instantaneously down to the water table — this means drape is redundant with newton. For what it's worth, this is not a permanent solution (we are considering more realistic options) but it is consistent with MODPATH 7.

We are digging in and hope to have a fix soon.

@aprovost-usgs please correct me if I got anything wrong here

[javgs-bd]

Hi @wpbonelli many thanks for looking into this. I just wanted to leave some comments to give more context on the particular case.

The example we built to report this is mimicking a more complex 3D tailings storage facility seepage model where only contaminant migration released into the initially dry more permeable layers are of interest. We are currently using MP7 because it is working as expected with newton, but we have to unrotate the whole grid for that. We are looking forward to migrate to PRT. I can't check today but I will do tomorrow. I'm almost sure that MP7 removes particles in dry cells when drape is on and the solution is using newton.

The drape action is critical because it prevents us to model many particles we don't need to. As the phreatic surface is raising, the model time and SP at which it reaches the initially dry layers is dependent on many the model parameters (namely, the hydraulic properties of different layers and tailings recharge rates), which are in turn, uncertain/adjustable parameters within an ensemble of models. We can't anticipate a precise release time, and therefore, removing particles released in dry layers makes a huge different in particle simulation times.

I wonder if leaving the released particles in dry layers with velocity = 0, as long as the head is below the cell bottom, could be an option in PRT. That would be neat (even better than removal via drape) because it would make the simulation require only one initial release. By releasing in many times or SPs, we are only trying to detect when the target layers became saturated and started contaminant movement.

[javgs-bd]

Hello @wpbonelli. Apologies for the delay on this. In this notebook we were checking that MP7 was eliminating particles released into dry layers even with newton activated. I hope this helps in the process of finding what's going on with PRT. As I commented previously, that v=0 for particles in released dry cells would be a life saver.

ex-prt-tr 3.zip

[wpbonelli]

@javgs-bd this model is very useful for debugging. Could we add it (with attribution) to the repo as a test case?