bug: flopy creates a column with "nan" values when assigning auxiliary variable in the EVT package

[giovannifi]

Describe the bug
I am building a flow+transpot model with Modflow 6 and I noticed a problem with the EVT package. I wanted to define 0 concentration to be extracted from the EVT. To do that I used the "auxiliary" variable option to be defined in the EVT package of the GWF model. By doing that I noticed that once the package is written there is, for each node where I define the evt, an additional and unwanted item at the end of each line with the value "nan". Is this the correct behaviour of the software?

To Reproduce
Steps to reproduce the behavior:

1. Please run the code below where I built a dummy GWF model (only the flow model). I used pycharm and python 3.9 to run the code below. In the code below I define only one node with evapotranspiration.

import flopy
import numpy as np

path_ws = r'F:\Working_Folder\testing\test_evt'
name = "test_evt"
h1 = 100
h2 = 90
Nlay = 10
N = 101
L = 400.0
H = 50.0
k = 1.0

sim = flopy.mf6.MFSimulation(sim_name=name, exe_name="mf6", version="mf6", sim_ws=path_ws)

tdis = flopy.mf6.ModflowTdis(sim, pname="tdis", time_units="DAYS", nper=1, perioddata=[(1.0, 1, 1.0)])

ims = flopy.mf6.ModflowIms(sim, pname="ims", complexity="SIMPLE")

model_nam_file = "{}.nam".format(name)

gwf = flopy.mf6.ModflowGwf(sim, modelname=name, model_nam_file=model_nam_file, )

bot = np.linspace(-H / Nlay, -H, Nlay)
delrow = delcol = L / (N - 1)
dis = flopy.mf6.ModflowGwfdis(gwf,
                              nlay=Nlay,
                              nrow=N,
                              ncol=N,
                              delr=delrow,
                              delc=delcol,
                              top=0.0,
                              botm=bot,
                              )

start = h1 * np.ones((Nlay, N, N))
ic = flopy.mf6.ModflowGwfic(gwf, pname="ic", strt=start)

npf = flopy.mf6.ModflowGwfnpf(gwf, icelltype=1, k=k, save_flows=True)

chd_rec = []
chd_rec.append(((0, int(N / 4), int(N / 4)), h2))
for layer in range(0, Nlay):
    for row_col in range(0, N):
        chd_rec.append(((layer, row_col, 0), h1))
        chd_rec.append(((layer, row_col, N - 1), h1))
        if row_col != 0 and row_col != N - 1:
            chd_rec.append(((layer, 0, row_col), h1))
            chd_rec.append(((layer, N - 1, row_col), h1))

chd = flopy.mf6.ModflowGwfchd(gwf,
                              maxbound=len(chd_rec),
                              stress_period_data=chd_rec,
                              save_flows=True,
                              )

iper = 0
ra = chd.stress_period_data.get_data(key=iper)

# Create the output control (`OC`) Package
headfile = "{}.hds".format(name)
head_filerecord = [headfile]
budgetfile = "{}.cbb".format(name)
budget_filerecord = [budgetfile]
saverecord = [("HEAD", "ALL"), ("BUDGET", "ALL")]
printrecord = [("HEAD", "LAST")]
oc = flopy.mf6.ModflowGwfoc(
    gwf,
    saverecord=saverecord,
    head_filerecord=head_filerecord,
    budget_filerecord=budget_filerecord,
    printrecord=printrecord,
)

rate = 0.001
ex_depth = 5.
concentration = 0.
evt = flopy.mf6.ModflowGwfevt(gwf,
                              maxbound=1,
                              nseg=1,
                              stress_period_data={0: [0, 0, 0, 0.0, rate, ex_depth, concentration]},
                              auxiliary='concentration'
                              )

sim.write_simulation()

2. open the generated EVT package and you will notice the "nan" values.

Expected behavior
I would expect the EVT to be created without the "nan" values.

Desktop (please complete the following information):

• OS: windows 10
• Browser: chrome
• python version: 3.9
• flopy version: 3.3.6

[langevin-usgs]

I think you are missing a surface elevation. You can see an example here.

[giovannifi]

In my case I set the surface elevation to 0. However, in my example I do not define PXDP and PETM because nseg = 1 and the manual says that "the PXDP and PETM arrays must be of size NSEG - 1". So my interpretation is that the length of the arrays for PXD and PETM should be 0 and therefore I do not have to define them. In your example you define NSEG = 1 and you have:

evtspd.append([(0, i, j), top[i, j], evtr, extdp, 1.0, 0.0])

I presume that the last term is the auxiliary variable (temperature) set to 0. However, I do not understand what 1 represents. Is it PXDP or PETM? and why is this defined?

On the basis of equation 6-33 of the documentation file called "tm6a55.pdf", if I have to declare something I should probably set:

PXDP =1 and PETM = 0

thanks for your support
Giovanni

[giovannifi]

Sorry for the multiple editing of my message, but I have just realised that in your example you defined SURF_RATE_SPECIFIED= True, this means that you had to define PETM0. So the "1" in your line should be for PETM0. However in my example I did not specify SURF_RATE_SPECIFIED= True, so I should not define anything.

[ccalvocm]

Greetings @giovannifi:

Please refer to the manual, there seems to be an issue in your stress_period_data. You are not missing the "surface" parameter, but might have mixed the order of the parameters in the stress_period_data dictionary which receives: cellid, surface, rate, depth, pxdp, petm, petm0, aux and boundname. You might get it if you rather write:

evt = flopy.mf6.ModflowGwfevt(gwf, maxbound=1, nseg=1, stress_period_data={0: [0, 0, rate, ex_depth,0, 0.0,0,concentration]}, auxiliary='concentration')

Kind regards

Tested on:

OS: windows 10
python version: 3.8.15
flopy version: 3.3.6

[giovannifi]

Hi @ccalvocm
thanks a lot for your reply. I still do not fully get it. I thought that pxdp, petm, petm0 were optional variables. The manual says that:

"the PXDP and PETM arrays must be of size NSEG - 1"

In my case nseg = 1, so PXDP and PETM should not be defined. In addition PETM0 is required only if the option "SURF_RATE_SPECIFIED " is set to True, which is not my case.

Indeed, when I define the EVT package without the transport component (without the auxiliary option) and I define the stress_period_data as follows:

stress_period_data = {0: [0, surf_elev, rate, ex_depth]

everything works fine. In addition, if it is like you said I should see at least three values with "nan" because I am missing the three parameters.

I still believe that the software is not behaving properly in this case.

Giovanni

[langevin-usgs]

HI @giovannifi, I took another look, and it does look like there is a problem with flopy. I've specified a list-based EVT record as shown below. In addition to the cellid, which is provided as a tuple, flopy insists that 5 additional items are entered, in addition to the cellid, even though surf_rate_specified=False.

evt_record = [(0, 0, 0), -1., -2., -3., -4., -5.]
evt = flopy.mf6.ModflowGwfevt(gwf,
                              maxbound=1,
                              nseg=1,
                              stress_period_data={0: [evt_record]},
                              auxiliary='concentration',
                              print_input=True,
                              surf_rate_specified=False,
                              )

The input file that is created looks like the following, which is consistent with the specified python input.

# File generated by Flopy version 3.3.7 on 05/12/2023 at 08:13:27.
BEGIN options
  auxiliary  concentration
  PRINT_INPUT
END options

BEGIN dimensions
  MAXBOUND  1
  NSEG  1
END dimensions

BEGIN period  1
  1 1 1      -1.00000000      -2.00000000      -3.00000000      -4.00000000      -5.00000000
END period  1

However, if we look at what MODFLOW 6 is doing, we see the problem. The following is from the gwf list file:

 ------------------------------------------------------------------------------------------------------
  EVT NO.    LAYER    ROW     COL       SURFACE         MAX. RATE        EXT. DEPTH     CONCENTRATION  
 ------------------------------------------------------------------------------------------------------
     1         1       1       1    -1.00000000      -2.00000000      -3.00000000      -4.00000000     
 ------------------------------------------------------------------------------------------------------

MODFLOW 6 reads 4 entries, in addition to the cellid.

We will take a look at this. @spaulins-usgs, do you think you could have a look?

[spaulins-usgs]

@langevin-usgs, confirmed that this is a problem. Currently the EVT dfn does not tie petm0 to surf_rate_specified. I will need to add a connection to these two variables. This has been typically done by setting one variable's "shape" based on another variable. In this case I will need to implement setting a variable's (petm0's) "shape" based on a boolean value (if surf_rate_specified is true, petm0 is of length 1, otherwise length 0).

[giovannifi]

Thanks for acknowledging the issue.