mf6: IRCH is repeated for subsequent stress periods in MODFLOW-6 but not in Flopy

[aleaf]

The documentation isn't very clear on this, but from talking to @langevin-usgs, IRCH only needs to be specified once, and then MODFLOW-6 will use that array for subsequent stress periods. However, in Flopy, specifying a single (2D) IRCH array results in a transient2d array that is populated with zeros after the initial stress period. The behavior is the same after loading a model with irch specified once.

import flopy
mf6 = flopy.mf6
sim_name = 'testsim'
model_name = 'testmodel'
out_dir = 'tmp'

tdis_name = '{}.tdis'.format(sim_name)
sim = mf6.MFSimulation(sim_name=sim_name,
                   version='mf6', exe_name=exe_name,
                   sim_ws=out_dir)
tdis_rc = [(6.0, 2, 1.0), (6.0, 3, 1.0), (6.0, 3, 1.0)]
tdis = mf6.ModflowTdis(sim, time_units='DAYS', nper=3,
                          perioddata=tdis_rc)

model = mf6.ModflowGwf(sim, modelname=model_name,
                         model_nam_file='{}.nam'.format(model_name))

dis = mf6.ModflowGwfdis(model, length_units='FEET', nlay=1,
                                     nrow=2, ncol=2, delr=500.0,
                                     delc=500.0,
                                     top=100.0, botm=50.0,
                                     filename='{}.dis'.format(model_name))
irch = [[0, 1], 
        [2, 3]]
rcha = mf6.ModflowGwfrcha(model, irch=irch, recharge={0: 1, 1: 1, 2: 1})
rcha.irch.array.sum(axis=(1, 2, 3))

Output:

array([6, 0, 0])

[spaulins-usgs]

If I understand correctly the desired behavior in the example above is for .array to return a transient2d array with dimensions (nper, nrow, ncol) but with all stress periods populated with the same array ([[0,1], [2,3]] repeated for every stress period). Is that correct? If so, I could implement a general rule for all transient2d arrays so that any stress period without data is populated with the the array from the closest stress period before it that has data. Would that general rule produce the desired behavior for all transient2d arrays in MF6?

[aleaf]

I would defer to @langevin-usgs. If IRCH can be specified by stress period, then it seems like it should be transient 2d. But then I think the behavior should be "what you see is what you get". If MODFLOW-6 is reusing IRCH for subsequent periods, then IRCH in Flopy should reflect that. Having the zeros in there makes it appear as though MODFLOW is assigning recharge to the top layer for those stress periods. The other issue I guess is with round-tripping. As it is now, if you were to simply copy IRCH from one recharge package to another, wouldn't the zeros then be applied to the second package? If the user originally specified IRCH as a 2D array, this wouldn't be the intended behavior. Thanks for fixing this.

[christianlangevin]

@aleaf, I took your example and ran it, and I'm not sure I see the behavior that you were indicating. If we write the rcha package from you script above, we get the following recharge input package

# File generated by Flopy version 3.3.3 on 12/10/2020 at 09:55:52.
BEGIN options
  READASARRAYS
END options

BEGIN period  1
  irch
    INTERNAL  FACTOR  1
    1  2
    3  4
  recharge
    CONSTANT       1.00000000
END period  1

BEGIN period  2
  recharge
    CONSTANT       1.00000000
END period  2

BEGIN period  3
  recharge
    CONSTANT       1.00000000
END period  3

It seems like this is the behavior that you want, right? IRCH is specified just once, and only for the first period.

I also see that if I try to access rcha.irch.array, then an error is returned:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-15-153adacc692a> in <module>
----> 1 rcha.irch.array

~/langevin/dev/flopy3fork.git/flopy/mf6/data/mfdata.py in array(self)
    267     def array(self):
    268         kwargs = {"array": True}
--> 269         return self.get_data(apply_mult=True, **kwargs)
    270 
    271     @property

~/langevin/dev/flopy3fork.git/flopy/mf6/data/mfdataarray.py in get_data(self, layer, apply_mult, **kwargs)
   1451                             output = data
   1452                         else:
-> 1453                             output = np.concatenate((output, data))
   1454                     return output
   1455                 else:

<__array_function__ internals> in concatenate(*args, **kwargs)

ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 4 dimension(s) and the array at index 1 has 2 dimension(s)

So @aleaf, I think what you are seeing is that the .array is telling you that the values are zero, but really they aren't. They are actually None.

Or maybe I'm missing something here? What should we expect to get back from rcha.irch.array (or MFTransientArray, more generally)? I think that is perhaps the issue here.

[spaulins-usgs]

@langevin-usgs I got the same result as @aleaf when I ran his script. I am not able to reproduce the error you are seeing. For "rcha.irch.array" I get what @aleaf is describing (see below). Flopy for MF6 transient array property currently functions similar to how the old flopy for MF2005 Transient2D array property functions, filling in values of an array of size [nper, nrow, ncol] and leaving any unfilled values equal to 0.

[[[[0 1]
[2 3]]]
[[[0 0]
[0 0]]]
[[[0 0]
[0 0]]]]

[jdhughes-dev]

@langevin-usgs @spaulins-usgs can this be closed?

[christianlangevin]

@spaulins-usgs, I believe you've started on this, or maybe even submitted a fix?

[spaulins-usgs]

@langevin-usgs We need to discuss the get_data method in the different data types to determine what the ideal behavior is for each data type. I sent you a spreadsheet outlining the current behavior which we need to discuss. I will resend and then let's set up a time to discuss.