Simulated water head higher than surface elevation

[Qing He]

Hi all,

I'm using modflow to simulate the groundwater head. The model runs successfully with converged water balance. However, when I checked the results I found the simulated water head shows to be higher than the land surface in some regions (but in other regions the results are normal and correspond well to the previous studies). 

I wonder what could be possible reasons for this? Any advice will be highly appreciated!

Best regards,

Qing

[Richard B. Winston]

1. The recharge rate may be too high.

2. The hydraulic conductivity is too low.

3. There could be boundaries that remove water but that were not included in the model.

[ashutosh singh]

Hi,

Just out of curiosity, can it be artesian conditions?

Not saying that you are modeling exactly as in my case. 

But I have seen it using simple modflow models that simulated groundwater levels that were higher than the surface elevation at places where open to air parking places existed.

To my surprise, I later found out that pumping wells were installed close by to keep the groundwater levels low because the phenomenon of groundwater uplift existed in the parking place.

These pumping wells were actually not there in my groundwater model.

//Ashutosh

--
This group was created in 2004 by Mr. C. P. Kumar, Former Scientist 'G', National Institute of Hydrology, Roorkee. Please visit his webpage at https://www.angelfire.com/nh/cpkumar/
---
You received this message because you are subscribed to the Google Groups "MODFLOW Users Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to modflow+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/modflow/d9bd7ec3-e722-4629-ab0b-3d25d757ae99%40mindspring.com.

[René M]

If it is allover the model space, I'd suggest you first look into your recharge (maybe too high) and K values (too low)!

A good appoach is to run the model with Zero Recharge and see if you still have high head. 

Otherwise, if the conditions occurs at part of the model, investigate if it is not due to artesian conditions.

//René

--

This group was created in 2004 by Mr. C. P. Kumar, Former Scientist 'G', National Institute of Hydrology, Roorkee. Please visit his webpage at https://www.angelfire.com/nh/cpkumar/
---
You received this message because you are subscribed to the Google Groups "MODFLOW Users Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to modflow+u...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/modflow/1213712d-39a1-4400-8602-72af95c55ce7n%40googlegroups.com.

[Qing He]

Thank you for all your replies. I feel it's not likely because of the missing artesians,  because there are vast areas showing water head higher than land surface. Almost over the majority of the eastern U.S. (see attached figure, non-shaded color indicate the irregular pixels). 

I checked the recharge input but it seems no problem. So I assume it might be caused by the low hydraulic conductivity. I'll do further check and get back to you when I found out the reason.

[Luther]

Could be that recharge is causing mounding. Lower hydraulic conductivity regions can make this issue worse. 

--

[Qing He]

I checked the conductivity data it seems no problem either (see the first figure below). In the areas where water head higher than land surface, the aquifer conductivity can be as high as 10 m/d. So I assume the problem of low conductivity can also be excluded. 

Could it be caused by the shallow aquifer thickness? The map I used has the thickness less than 20m in most areas (the second figure attached), but in some studies the thickness shows to be several hundred meters (see the third figure).

[Luther]

Yes, aquifer thickness is likely the problem, especially if the aquifer is much thicker in reality. The no-flow boundary at the model bottom and the less gradual change in specific discharge in the z direction causes the hydraulic head to artificially increase. 

To view this discussion on the web visit https://groups.google.com/d/msgid/modflow/1578ae7f-4165-46a0-a7db-04b4521d379cn%40googlegroups.com.

[René M]

Can you share you model configuration (e.g. how many layers, boundaries, etc.).

Another trick to curtail the high head is to apply a drain bc at say 1m below the surfave all over the model domain. But I feel the problem is still with your setup.

Have you tried the earlier suggestion to run the model with zero recharge? What was the result?

//

To view this discussion on the web visit https://groups.google.com/d/msgid/modflow/1578ae7f-4165-46a0-a7db-04b4521d379cn%40googlegroups.com.

[rbwinston]

Are you using consistent units throughout the model? For instance, if the grid cell size is in meters and the time step size is in seconds, the recharge and hydraulic conductivity should be in m/s. I am also concerned that in your image, Greenland and South America appear to be of similar sizes even though Greenland is much smaller than South America. If you aren't using an equal-area projection for your model, that could bias your results. 

To view this discussion on the web visit https://groups.google.com/d/msgid/modflow/1578ae7f-4165-46a0-a7db-04b4521d379cn%40googlegroups.com.

[Boyce, Scott E]

This is more of a hack until you get the properties dialed in, but you can always set up drains (DRN or DRT) at the land surface elevation. Any cells that have the head exceed the landsurface elevation will have the excess water removed. 

If you use MODFLOW-OWHM (https://code.usgs.gov/modflow/mf-owhm) you can use the Drain Return Flow package (DRT) and any water that is removed from the land surface via DRT can be passed to the Farm Process (FMP) as landscape runoff, which then ends up at SFR as runoff.

A common issue is that MODFLOW for convertible layers switches from Specific Yield to Specific Storage once the cell becomes fully saturated. If the top of the cell is the land surface, once it becomes fully saturated, the head will shoot upwards cause it switched form Sy to Ss. 

Another hack/fix is to have a large Ss for your upper most layer to keep the head from spiking above the land surface. Depending on your system it may even be a correct assumption cause the soil surface may not be a true confining barrier.

I'm not sure what version of MODFLOW you are using but packages like LPF and UPW take Specific Storage and multiply it by the aquifer thickness to get the Storativity. If you want to back door Specific Yield to be used even for the confine scenario, you take your Sy value and divide it by the aquifer thickness. Then when MODFLOW sets up the stortivity, it mimics Sy rather than Ss. So for example, say you want the confine situation to have a Sy = 0.3 and the top cell is 100 ft thick, then you would set Ss = Sy/thick = 0.003 and then when MODFLOW uses that value it multiplies the Ss by thick to get back the Sy.

But again, often when cells spike its because you have too much inflow somewhere, such as RCH, or a bad initial condition (such as lay 2 has too high of a head).

Hope that helps out,

Scott

---

From: mod...@googlegroups.com <mod...@googlegroups.com> on behalf of Qing He <heq...@g.ecc.u-tokyo.ac.jp>
Sent: Thursday, September 21, 2023 6:15 AM
To: MODFLOW Users Group <mod...@googlegroups.com>
Subject: [EXTERNAL] Re: [MODFLOW] Simulated water head higher than surface elevation

 

 

This email has been received from outside of DOI - Use caution before clicking on links, opening attachments, or responding.

To view this discussion on the web visit https://groups.google.com/d/msgid/modflow/1578ae7f-4165-46a0-a7db-04b4521d379cn%40googlegroups.com.

[Qing He]

Dear all,

I appreciate all your advice but let me ask some very general questions before conducting the sensitivity simulations as you suggested.

First, regarding the drainage boundary, I wonder how does the groundwater drains in MODFLOW? 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?

My next question is that could the initial conditions 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?

Then, related to the previous question, 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 a so called "spin-up" period for the model to reach the equilibrium state. But here as the model is running under the steady-state model, there is only one time step in total and no time for the model to spin-up. Should I change to the transient mode to see the results though?

I attached the model set-up 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")

[Qing He]

Dear all,

Thank you again for your very useful advice! I've now figured out what's the problem. It seems it's because the drainage and riverbed conductance i use are too low so it impedes the water removal from the system. I also figured out the answer about the other questions I've raised (initial condition and steady-state simulation). Thank you so so much again for your help!!