Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
javgs-bd opened this issue Sep 4, 2024 · 4 comments
Assignees
Labels
Milestone

Comments

@javgs-bd
Copy link

javgs-bd commented Sep 4, 2024

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.

  1. 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.

  1. 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
image
image
image

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
Copy link
Contributor

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 fix(prp): fix release time selection logic #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

wpbonelli added a commit that referenced this issue Sep 6, 2024
A conditional expression in the release timing logic was incorrect: we only want to enable the default release time (beginning of the first step of the first period) if no release times are specified via the RELEASETIMES block, the RELEASE_TIME_FREQUENCY option, or period block release time settings.

We also did not previously describe said default. Describe this in the PRP definition file and mention it in the release notes.

#2014 flushed this bug out, but this is not the main problem reported in that issue. More to come.
@javgs-bd
Copy link
Author

javgs-bd commented Sep 8, 2024

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
Copy link
Author

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

image

@wpbonelli wpbonelli added this to the 6.6.0 milestone Nov 20, 2024
@wpbonelli
Copy link
Contributor

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants