Problems with LAK and LKT packages: water and mass balance issues

[giovannifi]

I am trying to simulate a very simple model: I have a rectangular 1 layer (2d) model. The model is confined. I have a uniform flow set from left to right as steady state condition before applying the stresses that I will describe below.

After the steady state period, I set a drain and with the water mover package I route its volumes to a lake. The mover finally moves some water from the lake to an injection well. I need this particular set up, so do not suggest me to remove the lake and route the drain volumes directly to the injection well. The lake was defined only as a transfer pond with no interaction with the aquifer. I also defined the lake in a way that for each simulated period I have no storage (all the water in leaves the lake via an outlet that was set at the bottom of the lake).

The background concentration of the model is 20 g/L. However, I intentionally used the CNC package to define a constant concentration of 50 g/L at the location of the drain. My intention is to verify that the abstraction concentration from the drain is identical to the injection concentration at the well.

So, I created the model with the flopy code below. Please note that I simulated 1e6 days without the drain and injection as warm up period to reach steady state conditions (MF6 transport does not have steady state option).

If you run the code, please open the two files gwf_budget_lak.csv and gwt_budget_lak.csv generated after the simulatiuon. They are the files that record the water and mass balance at the lake respectively. In the first file you will recognise that the only non 0 values are:

INs: FROM-MVR_IN -> this is the volume of water from the Drain Package
OUTs: EXT-OUTFLOW_OUT, TO-MVR_OUT -> the sum of the two is the volume lost via the outlet. The second term is routed to the MAW package whilst the other is volume lost from the system

From the GWF point of view the above makes totally sense to me. However let's open the mass balance file of the lake (gwt_budget_lak.csv). In this case the non 0 values after the warm up period are:

INs: FROM-MVR_IN
OUTs: EXT-OUTFLOW_OUT, TO-MVR_OUT, STORAGE_OUT

I am not able to understand why I have the term STORAGE_OUT in this case and if you calculate the concentration of the water out (sum of EXT-OUTFLOW_OUT_GWT, TO-MVR_OUT_GWT / sum of EXT-OUTFLOW_OUT_GWF, TO-MVR_OUT_GWF) you will discover that it is variable with time and always lower than 50 g/l. However, the concentration of the water in is exactly 50 g/L as I was expecting.

Is anyone able to tell me why I get the term STORAGE_OUT in the LKT mass balance and why the concentration of the water out of the lake is not constant and equal to 50 g/L?

## Import Module Section
import numpy as np
import flopy

## Model Inputs
path_exe = r'C:\Program Files\Modflow\mf6.4.1\bin\mf6.exe'
path_model = r'F:\Working_Folder\testing\injection_opt\1_lake_transport_3'
modelname = 'lake'
run_mf = True

# Model domain and grid definition
Lx = 10000.
Ly = 10000.
ztop = 0.
zbot = -20.
nlay = 1
nrow = 100
ncol = 100
delr = Lx / ncol
delc = Ly / nrow
delv = (ztop - zbot) / nlay

# drain inputs
drn_r = 49
drn_c = 49
stage = -5.0
cond_drn = 100.0
idx_drn = (0, drn_r, drn_c)
conc_at_drains = 50.

# injection well inputs
wrow = 9
wcol = 19
idx_wel = (0, wrow, wcol)
head_limit = 30.5
qmax = 1000.
hk_skin = 1.e5

# temporal discretisation inputs
nper = 31
n_days = 10
n_ts = 100
# n_days_sp1 = 1
# n_ts_sp1 = 1
n_days_sp1 = 1e6
n_ts_sp1 = 100
# steady_state = True
steady_state = False
if not steady_state:
    sp_t = 0
else:
    sp_t = 1
        
# hydraulic properties of the model
kh = 10.0  # m/d
ss_value = 1.e-4

# constant heads
ch_east = 20.0
ch_west = 30.0
conc_west = 20.
conc_east = 0.

# initial head
h0 = 30.0

aux_name_conc = 'concentration'
conc0 = ch_east


## Model Build
# Create the Flopy simulation object
sim = flopy.mf6.MFSimulation(sim_name=modelname, exe_name=path_exe, sim_ws=path_model, continue_=True)

# Create the Flopy temporal discretization object
tdis_perioddata = [(n_days_sp1, n_ts_sp1, 1.0)]
for sp in range(nper-1):
    tdis_perioddata.append((n_days, n_ts, 1.0))
tdis = flopy.mf6.ModflowTdis(sim, pname='tdis', time_units='DAYS', nper=nper, perioddata=tdis_perioddata)

# Create the Flopy groundwater flow (gwf) model object
gwf_name = modelname + '_gwf'
gwf = flopy.mf6.ModflowGwf(sim, modelname=gwf_name, save_flows=True)

# discretization package
dis = flopy.mf6.ModflowGwfdis(gwf,
                              pname='dis',
                              nlay=nlay,
                              nrow=nrow,
                              ncol=ncol,
                              delr=delr,
                              delc=delc,
                              top=ztop,
                              botm=zbot)

# Create the initial conditions package
start = h0 * np.ones((nlay, nrow, ncol))
ic = flopy.mf6.ModflowGwfic(gwf, pname='ic', strt=start)

# Create the node property flow package
npf = flopy.mf6.ModflowGwfnpf(gwf, 
                              pname='npf', 
                              icelltype=0, 
                              k=kh, k33=kh, 
                              save_flows=True)
ss = flopy.mf6.ModflowGwfsto.ss.empty(gwf, 
                                      layered=True, 
                                      default_value=ss_value)
sto = flopy.mf6.ModflowGwfsto(gwf, 
                              pname='sto', 
                              save_flows=True, 
                              iconvert=1, ss=ss,
                              steady_state={0: steady_state}, 
                              transient={sp_t: True})

# Create the constant head package.
chd_rec = []
for r in range(nrow):
    chd_rec.append(((0, r, 0), ch_west, conc_west))
    chd_rec.append(((0, r, ncol - 1), ch_east, conc_east))
chd = flopy.mf6.modflow.mfgwfchd.ModflowGwfchd(gwf, 
                                               pname='chd', 
                                               maxbound=len(chd_rec),
                                               stress_period_data=chd_rec, 
                                               save_flows=True,
                                               auxiliary=aux_name_conc)

# drain package to simulate dewatering
drn_data = {}
drn_data[0] = [(idx_drn, stage, 0.)]
drn_data[1] = [(idx_drn, stage, cond_drn)]
n_drn = len(drn_data[0])

drn = flopy.mf6.ModflowGwfdrn(gwf, 
                              pname='drn-1', 
                              mover=True, 
                              save_flows=True, 
                              stress_period_data=drn_data,
                              maxbound=n_drn, 
                              filename=gwf_name + '_dr1.drn')

# drains at the well to control the max flow into the aquifer
drn_data2 = {}
drn_data2[0] = [(idx_wel, head_limit, 0.)]
drn_data2[1] = [(idx_wel, head_limit, 1.e5)]
n_drn2 = len(drn_data2[0])

drn2 = flopy.mf6.ModflowGwfdrn(gwf, 
                               pname='drn-2', 
                               save_flows=True, 
                               mover=True,
                               stress_period_data=drn_data2,
                               maxbound=n_drn2, 
                               filename=gwf_name + '_dr2.drn')

# LAK package to simulate transfer of water from drains to wells
idx_lak = (0, 20, 26)
bot_lake = 40.
str_lake = 40.
conc_lake = 0
lak_package_data = [[0, str_lake, 1, 'transfer_pond']]
lak_conndata = [(0, 0, idx_lak, 'VERTICAL', 0., 0, 0., 0., 0.)]
outlets = [(0, 0, -1, 'MANNING', bot_lake, 1e4, 1e-5, 0.5)]
# outlets = [(0, 0, -1, 'SPECIFIED', bot_lake, 10000, 1e-5, 0.1)]

lak_period_data = {}
lak_period_data[0] = [[0, 'STATUS', 'INACTIVE']]
lak_period_data[1] = [[0, 'STATUS', 'ACTIVE']]
# lak_period_data[1] = [[0, 'RATE', -qmax]]

obs_lak = [('lake_stage', 'stage', 1)]

lak = flopy.mf6.ModflowGwflak(gwf,
                              save_flows=True,
                              boundnames=True,
                              pname='lak-1',
                              nlakes=1,
                              noutlets=1,
                              mover=True,
                              packagedata=lak_package_data,
                              connectiondata=lak_conndata,
                              perioddata=lak_period_data,
                              outlets=outlets,
                              time_conversion=86400,
                              length_conversion=1.,
                              print_stage=True,
                              budgetcsv_filerecord='gwf_budget_lak.csv',
                              observations={'lake_stage.csv': obs_lak}
                              )

# MAW package
conc_well = 0.
well_packagedata = [(0, 0.01, zbot, h0, 'SPECIFIED', 1, conc_well)]
well_condata = [(0, 0, idx_wel, ztop, zbot, hk_skin, 0)]
well_perioddata = {0: [(0, 'STATUS', 'INACTIVE')]}

well_perioddata[1] = [(0, 'STATUS', 'ACTIVE')]


obs_wells = [('well1', 'head', 1)]
maw = flopy.mf6.ModflowGwfmaw(gwf,
                              pname='maw-1',
                              mover=True,
                              save_flows=True,
                              nmawwells=1,
                              no_well_storage=True,
                              packagedata=well_packagedata,
                              connectiondata=well_condata,
                              perioddata=well_perioddata,
                              auxiliary=[aux_name_conc],
                              observations={'obs_well1.csv': obs_wells}
                              )


# MVR package
mvr_perioddata = {}
mvr_perioddata[0] = []
mvr_perioddata[1] = [('drn-1', 0, 'lak-1', 0, 'EXCESS', 0.), 
                     ('lak-1', 0, 'maw-1', 0, 'UPTO', qmax)]
mvr_packages = [('drn-1',), ('lak-1',),  ('maw-1',)]

mvr = flopy.mf6.ModflowGwfmvr(gwf,
                              pname='mvr',
                              print_flows=True,
                              maxmvr=len(mvr_perioddata[1]),
                              maxpackages=len(mvr_packages),
                              packages=mvr_packages,
                              perioddata=mvr_perioddata,
                              budgetcsv_filerecord='gwf_budget_mvr.csv'
                              )

# output control (GWF)
oc = flopy.mf6.ModflowGwfoc(gwf, 
                            pname='oc', 
                            budget_filerecord='{}.cbb'.format(gwf_name),
                            head_filerecord='{}.hds'.format(gwf_name),
                            budgetcsv_filerecord = '{}.csv'.format(gwf_name),
                            saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')],
                           )

# IMS Package for the GWF
ims = flopy.mf6.ModflowIms(sim, pname='ims',
                           complexity='COMPLEX',
                           outer_dvclose=1e-10,
                           inner_dvclose=1e-10,
                          )
sim.register_ims_package(ims, [gwf.name])


# Transport model
gwt_name = modelname + '_gwt'
gwt = flopy.mf6.ModflowGwt(sim, modelname=gwt_name)

gwfgwt = flopy.mf6.ModflowGwfgwt(sim,
                                 exgtype='GWF6-GWT6',
                                 exgmnamea=gwf_name,
                                 exgmnameb=gwt_name)

# FMI package
# fmi = flopy.mf6.ModflowGwtfmi(gwt,
#                               flow_imbalance_correction=True)

# DISV package for the transport model
dis_t = flopy.mf6.ModflowGwtdis(gwt, 
                                nlay=nlay,
                                nrow=nrow,
                                ncol=ncol,
                                delr=delr,
                                delc=delc,
                                top=ztop,
                                botm=zbot)

# IC package (initial conditions concentration)
iconc = conc0 * np.ones((nlay, nrow, ncol))
iconc[0, drn_r, drn_c] = conc_at_drains
ic_t = flopy.mf6.ModflowGwtic(gwt, strt=iconc)

# ADV package
adv = flopy.mf6.ModflowGwtadv(gwt, scheme='tvd')

# DISP package
alh = 50.  # m
ath1 = alh / 10.
alv = ath1 / 10.
dsp = flopy.mf6.ModflowGwtdsp(gwt, alh=alh, ath1=ath1, alv=alv, diffc=0., xt3d_off=True)

# MST package
ne = 0.25
mst = flopy.mf6.ModflowGwtmst(gwt, porosity=ne)

# MWT package
mwt_packagedata = [0, conc0]
mwt_perioddata = {}
mwt_perioddata[0] = [(0, 'STATUS', 'INACTIVE')]
mwt_perioddata[1] = [(0, 'STATUS', 'ACTIVE')]
obs_wells_conc1 = [('well1_conc', aux_name_conc, 1)]

mwt1 = flopy.mf6.ModflowGwtmwt(gwt,
                               packagedata=mwt_packagedata,
                               mwtperioddata=mwt_perioddata,
                               flow_package_name='maw-1',
                               flow_package_auxiliary_name=aux_name_conc,
                               observations={'obs_conc_maw1.csv': obs_wells_conc1}
                               )

# SSM package
ssm_sources = [['chd', 'AUX', aux_name_conc]]
ssm = flopy.mf6.modflow.ModflowGwtssm(gwt,
                                      sources=ssm_sources,
                                      )

# MVT package
mvt = flopy.mf6.ModflowGwtmvt(gwt,
                              budgetcsv_filerecord='gwt_budget_mvt.csv'
                              )

# CNC package
cnc_period_data = {}
cnc_period_data[0] = [(idx_drn, conc_at_drains)]
cnc = flopy.mf6.ModflowGwtcnc(gwt,
                              maxbound=len(cnc_period_data[0]),
                              stress_period_data=cnc_period_data
                              )

# LKT package
lkt_packagedata=[(0, 0.)]
lkt_period_data = {}
lkt_period_data[0] = [[0, 'STATUS', 'INACTIVE']]
lkt_period_data[1] = [[0, 'STATUS', 'ACTIVE']]
obs_lkt = [('lake_conc', 'CONCENTRATION', 1)]

lkt = flopy.mf6.ModflowGwtlkt(gwt,
                              flow_package_name='lak-1',
                              packagedata=lkt_packagedata,
                              lakeperioddata=lkt_period_data,
                              observations={'lake_conc.csv': obs_lkt},
                              budgetcsv_filerecord='gwt_budget_lak.csv',
                             )
    
# OC package (GWT)
concfile = '{}.ucn'.format(gwt_name)
conc_filerecord = [concfile]
saverecord_t = {0: [('CONCENTRATION', "ALL")]}
printrecord_t = [('CONCENTRATION', "LAST")]
oc_t = flopy.mf6.ModflowGwtoc(gwt,
                              saverecord=saverecord_t,
                              concentration_filerecord=conc_filerecord,
                              printrecord=printrecord_t,
                              budgetcsv_filerecord='{}.csv'.format(gwt_name),
                              )


# IMS Package for GWT
imst = flopy.mf6.ModflowIms(sim,
                            filename=modelname + '_gwt.ims',
                            csv_outer_output_filerecord='ims_outer_output_filerecord_t.csv',
                            csv_inner_output_filerecord='ims_inner_output_filerecord_t.csv',
                            complexity='COMPLEX',
                            outer_dvclose=0.0001,
                            inner_dvclose=0.0001,
                            relaxation_factor=0.97,
                            outer_maximum=100,
                            inner_maximum=300,
                            inner_rclose=10000.
                            )
sim.register_ims_package(imst, [gwt.name])


## Write packages and run the model
sim.write_simulation()

# Run the simulation
if run_mf:
    success, buff = sim.run_simulation(silent=True)
    print('Success is: ' + str(success))
else:
    success = True

[langevin-usgs]

You've set the initial concentration of the lake to be zero in the lkt_packagedata. The LKT Package looks at the solute mass inflow (from the drain) and the outflows (ext-outflow, to-mvr, and storage) and recalculates the lake concentration. The LKT storage term you are seeing is not a flow storage term, but rather the amount of mass in the lake that has increased or decreased during the time step. Any chance this explains what you are seeing? Maybe if you set the initial concentration of the lake to be 50, the model will behave as expected?

[giovannifi]

@langevin-usgs thanks a lot for your prompt reply, but I am not able to fully understand your explanation and I believe that your explanation is not consistent with the logic of the flow component, i.e. why do not we have a "storage" component in the GWF balance of the lake? I. e. why do not we have a component that accounts for "the amount of [water] in the lake that has increased or decreased during the time step"? why is the logic different for the mass?

In reality, I believe that this does not make any sense, the lake does not store any water and there could not be any increase or decrease of mass.

Indeed, I intentionally set the initial concentration of the lak to 0. But I also set the initial stage of the lake to 0 as well, so no volume of water is present in the lake at the beginning. The lake is also set in a way that any ingress of water leaves the lake via the outlet without allowing any storage. Therefore, one would expect that the concentration of the water in should be the same as the concentration of the water out. Still do not fully understand the "storage" term in the mass balance.

If I set the initial concentration of the lake to 50, I would completely defeat the purpose to know the concentration of the abstracted waters from the drains. In this particular case I control it with the CNC package and it is easy to guess the concentration, but in my real case scenario this does not happen and I was simply testing this synthetic case to understand if the mover was able to route the right concentration to the injection borefield. If this works, this would avoid me to intervene during the simulation with the API. Essentially I am looking for a method that is able to route the right concentration from the drain to the injection well. This works great for the flow component, but, to me, it does not seem to work for the mass.

As mentioned in my previous message, I am forced to use the lake as transfer pond, because in my real case scenario I have to route the abstracted waters from the drains to an injection borefield made of multiple bores.

Does it make sense?

thanks again for your help

[giovannifi]

@langevin-usgs sorry for my message above, but I found my error. I actually set the initial stage of the lake to be 40 rather than 0. So my lake is starting with unwanted initial volume of water. This explains the storage term now and it makes totally sense. I re-run the script by changing the variable str_lake = 0. I also noticed that I had to set the variable INVERT for the outlet to be 0 as well in order to have the outlet at the bottom of the lake. I thought this should have been the absolute elevation and not the stage. With these modifications I get the concentration I wanted for the water out (from the lake).