Under relaxation: pumping dewatering cell

[ougx]

I was testing the under relaxation option. Below is the code that I created two models: one is a one-column model and the other is a two-column model. Both model with one pumping well that will dewater the top cell. The newton and under relaxation are on. While everything else is the same, the two-column model cannot converge. Seems to me under relaxation did not work for the second model.

import flopy
import os

nlay = 2
nrow = 1
ncol = 1
thick = 1
dx = 1
dy = 1

modelname = 'L2C1'

top = nlay*thick
botm = [[[thick * (nlay - i -1)]] for i in range(nlay)]

def write_col1(sim_ws='.'):
    os.makedirs(sim_ws, exist_ok=True)
    sim = flopy.mf6.MFSimulation(sim_ws=sim_ws)
    gwf = flopy.mf6.MFModel(sim, model_type='gwf6', modelname=modelname, )
    gwf.name_file.newtonoptions = "NEWTON UNDER_RELAXATION"
    tdis = flopy.mf6.ModflowTdis(sim, perioddata=[[1, 1, 1.]])
    dis = flopy.mf6.ModflowGwfdis(gwf, nlay=nlay, nrow=nrow, ncol=ncol, delr=dx, delc=dy, top=top, botm=botm, idomain=[1,0])
    npf = flopy.mf6.ModflowGwfnpf(gwf, icelltype=1, k=1, k33=0.1)
    sto = flopy.mf6.ModflowGwfsto(gwf, iconvert =1, ss=1e-5, sy=0.1)
    wel = flopy.mf6.ModflowGwfwel(gwf, stress_period_data=[[0,0,0,-0.101]])
    ims = flopy.mf6.ModflowIms(sim, print_option='summary', outer_maximum=100, inner_maximum=100, linear_acceleration="BICGSTAB")
    ic = flopy.mf6.ModflowGwfic(gwf, strt=top,)
    oc = flopy.mf6.ModflowGwfoc(gwf, head_filerecord=modelname+'.hds', saverecord=['HEAD', 'all'], printrecord=['HEAD', 'all'])
    sim.register_ims_package(ims, modelname)
    sim.write_simulation()
write_col1('col1')


ncol = 2
botm = [[[thick * (nlay - i -1),]*ncol]*nrow for i in range(nlay)]
#%% 
def write_col2(sim_ws='.'):
    os.makedirs(sim_ws, exist_ok=True)
    sim = flopy.mf6.MFSimulation(sim_ws=sim_ws)
    gwf = flopy.mf6.MFModel(sim, model_type='gwf6', modelname=modelname, )
    gwf.name_file.newtonoptions = "NEWTON UNDER_RELAXATION"
    tdis = flopy.mf6.ModflowTdis(sim, perioddata=[[1, 1, 1.]])
    dis = flopy.mf6.ModflowGwfdis(gwf, nlay=nlay, nrow=nrow, ncol=ncol, delr=dx, delc=dy, top=top, botm=botm, idomain=[[[1,1]],[[0,1]]])
    npf = flopy.mf6.ModflowGwfnpf(gwf, icelltype=1, k=1e-4, k33=0.1)
    sto = flopy.mf6.ModflowGwfsto(gwf, iconvert =1, ss=1e-5, sy=0.1)
    wel = flopy.mf6.ModflowGwfwel(gwf, stress_period_data=[[0,0,0,-0.101]])
    ims = flopy.mf6.ModflowIms(sim, print_option='summary', outer_maximum=100, inner_maximum=100, linear_acceleration="BICGSTAB")
    ic = flopy.mf6.ModflowGwfic(gwf, strt=top,)
    oc = flopy.mf6.ModflowGwfoc(gwf, head_filerecord=modelname+'.hds', saverecord=['HEAD', 'all'], printrecord=['HEAD', 'all'])
    sim.register_ims_package(ims, modelname)
    sim.write_simulation()
write_col2('col2')

[vivekbedekar]

The problem is designed such that the well pumps a quantity just more than the cell volume, causing the water level to drop below the bottom of layer 1. You are correct that one of the models does not converge. It looks like underrelaxation is not active in both the models. Why do you suspect that underrelaxation is the problem? Maybe you have a version in which underrelaxation was active and that caused the issue? It appears that the non-convergence is simply because of the non-linearity of the two-cell problem. Low K in layer 1 between the well cell and the neighboring cell. If the K is set to any higher value or a lower value, the model converges.