Issues in simulating a confined, infinitely extended aquifer.

[Guglio89]

Hi everyone, I hope you can support me. I am trying to simulate a confined, homogeneous, infinitely extended aquifer with a pumping well at its center under steady-state conditions. The pumping rate is constant.
The problem is as follows: when running simulations with Lx=Ly=200, nrow=ncol=200, delr=delc=1 under these conditions, I obtain a head.min() value of -3.0045915. However, when setting Lx=Ly=800, nrow=ncol=800, delr=delc=1, the value of head.min() is -3.6683185. I expect that, given the stationary conditions and constant pumping, and considering an infinitely extended aquifer, at a certain distance from the well there should be no impact from pumping, and that as the number of cells and the size of the aquifer increase, the values of head.min() should not change, provided that the values of delr and delc remain constant. The code is the following:
`import os
import numpy as np
import matplotlib.pyplot as plt
import flopy
import flopy.utils.binaryfile as bf

modelname = "Homogeneous"
mf = flopy.modflow.Modflow(
modelname, exe_name=r'C:\Users\gugli\Desktop\MODFLOW2005\MF2005.1_12\bin\mf2005.exe',
model_ws=r'C:\Users\gugli\Desktop\MODFLOW_analysis\Test')

#DIS PACKAGE
Lx =200.0
Ly = 200.0
ztop = 0.0
zbot = -20.0
nlay = 1
nrow = 200
ncol = 200
delr = Lx / ncol
delc = Ly / nrow
delv = (ztop - zbot) / nlay
botm = np.linspace(ztop, zbot, nlay + 1)
lenuni = 2 #Length units = meters

nper = 1
itmuni = 4 #time units = Days
perlen = [1]
nstp = [1]
steady = [True]

dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc, top=ztop, botm=botm[1:], nper=nper, perlen=perlen, nstp=nstp, steady=steady, itmuni=itmuni, lenuni = lenuni)

#BAS PACKAGE
ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)
ibound[:, :, 0] = -1
ibound[:, :, -1] = -1
ibound[:, 0, :] = -1
ibound[:, -1, :] = -1

strt = np.ones((nlay, nrow, ncol), dtype=np.float32) * 0

bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)

#GROUNDWATER FLOW PACKAGE: LAYER PROPERTY FLOW PACKAGE

Define hydraulic conductivity and specific storage

hk = 50 # hydraulic conductivity (m/day)
ss = 0.0001 # specific storage (1/m)
lpf = flopy.modflow.ModflowLpf(mf, laytyp=0, hk=hk, ss=ss, ipakcb=53)

#TRANSIENT WELL PACKAGE
pumping_rate = -3000.0 #m^3/day
wel_rec = [0, nrow / 2 - 1, ncol / 2 - 1, pumping_rate]
wel = flopy.modflow.ModflowWel(mf, stress_period_data=wel_rec)

#OUTPUT CONTROL
spd = {(0, 0): ["print head", "print budget", "save head", "save budget"]}
oc = flopy.modflow.ModflowOc(mf, stress_period_data=spd, compact=False)

#SOLVER:PRECONDITIONED CONJUGATE-GRADIENT PACKAGE

Define solver parameters

nper = 1 # number of stress periods
steady = [True] # steady-state stress period
flux_convergence = 0.001 # flux convergence criteria
max_iter = 100 # maximum number of iterations
head_tolerance = 0.0001 # head change tolerance criteria

pcg = flopy.modflow.ModflowPcg(mf, hclose=flux_convergence, rclose=head_tolerance, mxiter=max_iter)

#INPUT WRITE
mf.write_input()

#RUN
success, mfoutput = mf.run_model(silent=True, pause=False)
assert success, "MODFLOW did not terminate normally."

"POST-PROCESSING THE RESULTS"
#GET HEADS
headobj = bf.HeadFile(modelname + ".hds")
times = headobj.get_times()
head = headobj.get_data(totim=times[-1])`

[langevin-usgs]

If you are simulating steady state conditions with a pumping well, then all of the water withdrawn from the well must come from a boundary. Thus, for steady state conditions, you will see the effect of the well extend to the boundary, no matter how far away it is located from the well. In your case, you are surrounding your model with constant heads (with ibound=-1), and thus, the inflow into the model from the constant heads should equal the pumping rate. At your constant head boundaries, the drawdown is zero, by definition. I would expect the drawdown in the well to decrease as you move your boundaries farther away from the well.

I think part of the issue is that there isn't an answer to your problem. The steady-state solution for pumping from an infinitely extended aquifer is undefined. Note that a transient analysis for this problem will give very different results as water can also come from storage.

[Guglio89]

Thank you very much for your help. I just have one last doubt. If I keep the number of ncol and nrow cells constant, regardless of the values of Lx and Ly, and consequently of delr and delc, the values of head.min() do not change. For example, if for Lx=Ly=800, nrow=ncol=800, delr=delc=1, the value of head.min() is -3.6683185, then by setting Lx=Ly=8000, nrow=ncol=800, delr=delc=10, the value of head.min() still remains -3.6683185. I am not sure if my interpretation is correct, but it seems that it is not the size of the aquifer or the cell dimensions, but rather their number, that influences the result. At least for the case of a homogeneous aquifer. PS: I had thought that a solution to the problem could be to set ibound=1 everywhere, so that modflow could calculate the heads throughout the domain. However, the simulation does not terminate correctly.