A transient model with ~300 linesinks not matching 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.