General questions about how the steady-state simulation works in modflow

[Qing-sudo]

Hi all,

I'm using FloPy to do a 40year steady-state simulation but got the result in which the simulated groundwater head is much higher than the land surface over the lowlands (mountainous areas are okay). I got some feedbacks that there might be some bugs in my model set-up, but since I'm still very new to this model, I want to ask some general question about how the steady-state simulation works in modflow.

The most likely reason for this water-head-higher-than-land-surface problem could be the absence of water-removal terms in my model. I already included a drain package in my model but I wonder does the drain package can only simulate the groundwater flux to river? Theoretically the fate of groundwater flux should be ocean or the terminal rivers/lakes. I considered rivers and wells in my model, but the water head over ocean are set as constant 0m throughout the simulation. Could this cause limited amount of water-removal, and eventually, water head higher than land surface?

Another possible reason could be that the initial conditions can influence the water head. Currently I set the initial water head to 0m at all grid cells. Does this mean the water fluxes will accumulate based on this level, and therefore the final water head will be higher than the land surface over the lowlands?

My next question is not directly related to the problem, but since the simulation was conducted under the steady-state mode and there is only one time step in total, I wonder how does MODFLOW conduct the equilibrium simulation? Because if the initial water head is 0 all over the grids but the final water head is heterogenous in my study domain, there should be several time steps conducted as the so called "spin-up" period for the model to reach the equilibrium state. Should I change to the transient mode to see the results though?

I also attached the FloPy code below:

#========================build model============================#
prefix  = '/Volumes/heqing_drive/projects_backup/mf_glb_build/'
modelname = 'GLM40yr_stest2'   #maximum 16 characters
modelws   = prefix + modelname   #subdir where the output is going to be stored

#========================Read in geo data=======================#
nn = geoinput_var[0].shape
#define model parameters
#geological parameters of the study area: horizontal
ncol = nn[1]      #number of columns
nrow = nn[0]      #number of rows

delc = 111*(1/12)*1000   #meridional distance of each grid, 5min ~ 92.5km
delr = 111*(1/12)*1000   #zonal distance of each grid, 5min

#geological parameters of the study area: vertical (from GLHYMPS and HydroSHEDS)
nlay = 1       #total aquifer layers

top  = geoinput_var[0]
dtb  = geoinput_var[1]
botm = top - dtb

#aquifer property parameters (from GLHYMPS)
k_hor = geoinput_var[2]
k_ver = k_hor

#river geometry parameters
rbot = geoinput_var[4]
hriv = geoinput_var[5]
rb_cond = geoinput_var[6]

#recharge
rch_h08 = geoinput_var[7]/1000   #from mm/d to m/d

#pumping
pump_h08 = geoinput_var[8]

#==========================create a simulation============================#
sim = fp.mf6.MFSimulation(sim_name = modelname,
                         version  = 'mf6',
                         exe_name = '/Users/qinghe/modflow6/bin/mf6',   
                         sim_ws   = modelws)

#===========================time discretization===========================#
tdis_par = [(40*365.0, 1, 1.0)]   #simlation period length, time step, scale factor
tdis = fp.mf6.ModflowTdis(sim,
                         units = 'Days',
                         perioddata = tdis_par)

#===============interative model solution package (mandatory for all simulation cases)==========#
ims = fp.mf6.ModflowIms(sim,
                        complexity = 'SIMPLE',                        
                        outer_dvclose= 1e-3, # defining the criteria for maximum head change (outer iteration); default value is 1e-4
                        inner_dvclose= 1e-3, # defining the criteria for maximum head change (inner iteration); default value is 1e-5
                        rcloserecord = 0.1,  #In-Out criteria, default = 0.1
                        inner_maximum = 500, # defining the maximum number of inner (linear) iterations; default value is 50;  For nonlinear problems, I              NNER_MAXIMUM usually ranges from 60 to 600; a value of 100 will be sufficient for most linear problems.
                        outer_maximum = 100,
                        print_option = 'ALL')

#=============================groundwater flow package========================#
gwf = fp.mf6.ModflowGwf(sim,
                       modelname = modelname,
                       model_nam_file = f"{modelname}.name",
                       save_flows = True)

#=============================spatial discretization package======================#
ddis = fp.mf6.ModflowGwfdis(gwf,
                           length_units = "METERS",
                           nlay = nlay,
                           nrow = nrow,
                           ncol = ncol,
                           delr = delr,
                           delc = delc,
                           top = top,
                           botm = botm)

#=============================node property packag3===============================#
npf = fp.mf6.ModflowGwfnpf(gwf,
                          save_specific_discharge = True,
                          icellType = 1,
                          k = k_hor,     #hydro cond in x-axis direction
                          k22 = k_hor,   #hydro cond in y-axis direction
                          k33 = k_ver)   #hydro cond in z-axis direction

#=================initial water head for all grids (mendatory for all simulation cases)===============#
init_head = 0
ic = fp.mf6.ModflowGwfic(gwf,
                        strt = init_head)

#====================add constant head to sea areas (boundary conditions)====================#
#0m constant water head for sea (boundary condition)
constant_head = []
for i in range(len(sea_row)):
    constant_head.append([(0, sea_row[i], sea_col[i]), 0])   
       
print("Adding ", len(constant_head), "Constant heads for boundary condition")

#add constant head to sea areas (boundary conditions)
chd = fp.mf6.ModflowGwfchd(gwf,
                          stress_period_data = {0:constant_head})

#==============================for river package===============================#
riv_head = []
for i in range(len(riv_row)):
   
    ihriv  = hriv[riv_row[i],riv_col[i]]
    ircond = rb_cond[riv_row[i],riv_col[i]]
    irbot  = rbot[riv_row[i],riv_col[i]]
    riv_head.append([0, riv_row[i], riv_col[i], ihriv, ircond, irbot])
   
riv = fp.mf6.ModflowGwfriv(gwf,
                          stress_period_data = {0:riv_head})

print("Adding ", len(riv_head), "River heads for boundary condition")

#=============================== for drainage ===============================#
drn_elv = []
for i in range(len(drn_row)):
   
    ielv  = top[drn_row[i],drn_col[i]]
    icond = rb_cond[drn_row[i],drn_col[i]]
    drn_elv.append([0,drn_row[i],drn_col[i], ielv, icond])
   
drn = fp.mf6.ModflowGwfdrn(gwf,
                           stress_period_data = {0:drn_elv})
print("Adding ", len(drn_elv), "Drainage elevation for boundary condition")

#==============================for recharge package==========================#
rch = fp.mf6.ModflowGwfrcha(gwf,
                            recharge = rch_h08,
                            pname = 'rch')

#==============================for pumping package=============================#
pump_rate = []
for i in range(len(pump_row)):     #extract row and col for each well
    irate = -pump_h08[pump_row[i],pump_col[i]]   #negative indicates extraction
    pump_rate.append([(0, pump_row[i], pump_col[i]), irate])   #assign discharge data for each well; 3-D location of each well is (ilay,irow,icol)

wel = fp.mf6.ModflowGwfwel(gwf,
                           stress_period_data = pump_rate)
print("Adding ", len(pump_rate), "Pumping rate for boundary condition")

[RyanConway91]

If the heads are above the ground surface, you need something to remove water. A common approach is to use DRNs at the elevation of the ground surface. These represent water that would be lost to overland flow/surface water. Ideally the flux that ends up in these should make sense based on what you know about the regional water balance.

When solving a steady-state simulation, the initial heads should not matter if you're model is converging well. In a SS simulation, all the initial heads do is provide an initial condition to start trying to iteratively solve the GWF equation. It is common to set them to ground surface, as this is easy and often leads to decently quick convergence.