How to make a transient model with ~300 linesinks match MF6

[ougx]

Hi @dbrakenhoff and @mbakker7:

I followed the meandering river example to test its capability for stream depletion. Below is my script and data is attached riv1.csv.

The simulation is pumping near a stream for 50 yrs with a constant rate. The parameters are listed below. grid.csv includes the simulated drawdown with MF6.

T = 1000 ft^2/d
S = 0.1
Q = 20000 ft^3/d

import pandas as pd
import timflow.transient as tft

grid = pd.read_csv("grid.csv")
riv = pd.read_csv("riv1.csv")

# Create model
yr50 = 365.25*50
ml = tft.ModelMaq(
    kaq=10, z=[100, 0], Saq=[0.1], phreatictop=False, tmin=yr50, tmax=yr50 # 50 years
)

# Add a river with a fixed water level
river1 = tft.RiverString(ml, xy=riv.values, tsandh='fixed')

# Add a well
xw =2485446.31
yw = 441303.66
well = tft.Well(ml, xw, yw, rw=0.5, tsandQ=[(0, 20000)])

# Solve model
ml.solve()

# Compare with MF6
heads = []
fails = []
for i, r in grid.iterrows():
    try:
        grid.loc[i, "tim1"] = ml.head(r.x, r.y, yr50)
    except:
        pass
    

ax = grid.set_index("mf6").tim1.plot(linewidth=0, marker=".")
ax.grid()

timflow underpredicted the drawdown significantly compared with MF6. Did I do something wrong here? Appreciate your help.

BTW, I have signed up for the short coarse at MODFLOW & More. See you soon.

[mbakker7]

Thanks for your question, @ougx
You specified phreatictop=False and a specifc storage Saq=0.1. So what timflow does is, it will compute the storage coefficient as Saq multiplied by the aquifer thickness. And I don't think that is what you want (for your case it will compute a storage coefficient of 10!) So what you should do is phreatictop=True.

In the latest version of timflow, the phreatictop keyword is deprecated and you have to specify topboundary='phreatic'. That way timflow knows that the first value of Saq is phreatic storage and it won't multiply it with the aquifer thickness to obtain the storage coefficient of the aquifer.

The other thing I noticed is that your tmin and tmax are the same. tmin is the minimum time at which you can compute a head after a change in boundary condition. In your case the boundary condition changes at t=0 (when the well turns on). So in theory it should go ok. timflow arranges time by log cycles so tmin=10_000 and tmax=100_000 would be better (so that your time of interest is somewhere in that interval).

Does that resolve your differences?

[mbakker7]

I am going to move this to Discussions.

[mbakker7]

If you could share the MF6 grid results so that I can contour the results (or a short Python script that reads the grid.csv file and contours the values) that would make for a nice comparison.

[ougx]

Thanks Mark! They are well agreed after I changed Saq to 0.001. Attached is the MF6 files mf.zip. It use DISV grid so you could read the grid from the DISV file using flopy (I am sure you know this).

Noted that I have increased the complexity by adding another segment riv2.csv. Below is the updated script.

I used a large tmin because it takes very long to solve when tmin is small. Is it possible to extract the cumulative river discharge?

import pandas as pd
import timflow.transient as tft

grid = pd.read_csv("grid.csv")
riv1 = pd.read_csv("riv1.csv")
riv2 = pd.read_csv("riv2.csv")

# Create model
yr50 = 365.25*50
ml = tft.ModelMaq(
    kaq=10, z=[100, 0], Saq=[1e-3], phreatictop=False, tmin=10_000, tmax=100_000 # 50 years
)

# Add a river with a fixed water level
river1 = tft.RiverString(ml, xy=riv1.values, tsandh='fixed')
river2 = tft.RiverString(ml, xy=riv2.values, tsandh='fixed')

# Add a well
xw =2485446.31
yw = 441303.66
well = tft.Well(ml, xw, yw, rw=0.5, tsandQ=[(0, 20000)])

# Solve model
ml.solve()

# extract head
for i, r in grid.iterrows():
    try:
        grid.loc[i, "tim1"] = ml.head(r.x, r.y, yr50)
    except:
        pass


ax = grid.set_index("mf6").tim1.plot(linewidth=0, marker=".")
lims = [-31, 1.05]
ax.plot(lims, lims, "k--", lw=1, alpha=0.6, label="1:1 line")
ax.legend()
ax.grid(alpha=0.3, lw=0.5)

[mbakker7]

Thanks. I will try to make a comparison plot.

You are right that with a small tmin it will take longer to compute as a separate solution is obtained for every log cycle of time (but a tmin=10_000 will be just as quick for your case as it is still the same log cycle).

Regarding the river discharge, as you called your river string river1, you can use the river1.discharge(t) function, where t is time (and can be an array as long as the t values are between your tmin and tmax). The .discharge function gives the discharge at that time, not the cumulative discharge. There is not function for that, so you would have to compute the .discharge function every year or so for 50 years and then sum them up.

One more thing: You can speed up your model by using fewer River segments (you don't need 300). You can use a function from numpy or scipy to interpolate to say 50 elements and get pretty much the same answer. Also, if you want to compute values on a regular grid, you can use the headgrid or plots.contour function, which now have a parallel keyword argument. When that is set to True it will use all the available cores on your machine to compute a grid, which may speed up things.

[ougx]

Thanks, Mark. That will be another useful test case: evaluating how much we can simplify the stream geometry representation in Timflow without compromising performance. At the moment, I am benchmarking Timflow for stream depletion, so both Timflow and MODFLOW are using the same setup for comparison.