Coastal aquifers and the interaction among brine / fresh groundwater need to be evaluated with a modeling code that can deal with variable density flow. For more than 18 years, SEAWAT was the prefered (or only) open source solution implemented in Modflow 2000 with some limitations* on its use with Flopy. Now in Modflow 6 the concept of simulation involves flow and transport modeling together with exchange among them. We have done a tutorial with a simple case of variable density flow from a saline lake into an acuifer. The transient model has a duration of 50 days where the saline water "intrudes" the aquifer at the bottom part of the lake.

* When we talk about limitations, we talk about the gap on the capacities from the general public to deal with intermediate/complex Python code.

Tutorial

Python code

The code is divided in two parts: model simulation and output analysis.

Model simulation

This is the Python code for the model simulation.

#import required packages
import os
import flopy
import numpy as np
import matplotlib.pyplot as plt

#horizontal discretization
lx = 7.
lz = 4.
nrow = 1
ncol = 7
delc = 1.
delr = lx / ncol


#vertical discretization
nlay = 4
top = 4.
botm = [3., 2., 1., 0.]
delz = lz / nlay

#temporal discretization
nper = 1
perlen = 10.0
nstp = 50
tsmult = 1.

MF6 Simulation

ws = os.path.abspath('../Model/')
simName = 'buyExample'
sim = flopy.mf6.MFSimulation(sim_name=simName, version='mf6',
                             exe_name='../Exe/mf6',sim_ws=ws)

tdis = flopy.mf6.ModflowTdis(sim, time_units='DAYS',
                                 nper=nper, perioddata=[[perlen,nstp,tsmult]])

Flow model (gwf)

name = 'buyLake'
fModName = 'buyFlowModel'
tModName = 'buyTrspModel'
gwf = flopy.mf6.ModflowGwf(sim, modelname=fModName, newtonoptions=True)

nouter, ninner = 700, 300
hclose, rclose, relax = 1e-8, 1e-6, 0.97
imsgwf = flopy.mf6.ModflowIms(sim, print_option='ALL',
                                  outer_dvclose=hclose,
                                  outer_maximum=nouter,
                                  under_relaxation='NONE',
                                  inner_maximum=ninner,
                                  inner_dvclose=hclose, rcloserecord=rclose,
                                  linear_acceleration='BICGSTAB',
                                  scaling_method='NONE',
                                  reordering_method='NONE',
                                  relaxation_factor=relax,
                                  filename='{}.ims'.format(fModName))

idomain = np.full((nlay, nrow, ncol), 1)
idomain[0, 0, 1:6] = 0
idomain[1, 0, 2:5] = 0
idomain[2, 0, 3:4] = 0
dis = flopy.mf6.ModflowGwfdis(gwf, nlay=nlay, nrow=nrow, ncol=ncol,
                              delr=delr, delc=delc,
                              top=top, botm=botm, idomain=idomain)

#plot aerial grid
fig, ax = plt.subplots(figsize=(14,6))
mapview = flopy.plot.PlotMapView(model=gwf)
linecollection = mapview.plot_grid(alpha=0.25)

#plot cross section with active zone
fig, ax = plt.subplots(figsize=(14,6))
crossview = flopy.plot.PlotCrossSection(model=gwf, line={'row':0})
crossview.plot_grid(alpha=0.25)
crossview.plot_ibound()

<matplotlib.collections.PatchCollection at 0x1c2b67e74c0>

strt = np.zeros((nlay, nrow, ncol), dtype=np.float)
strt[0, 0, :] = 3.5
strt[1, 0, :] = 3.0
strt[1, 0, 1:6] = 2.5
strt[2, 0, :] = 2.0
strt[3, 0, :] = 1.0
ic = flopy.mf6.ModflowGwfic(gwf, strt=strt)

#plot cross section with initial heads
fig, ax = plt.subplots(figsize=(18,6))
crossview = flopy.plot.PlotCrossSection(model=gwf, line={'row':0})
crossview.plot_grid(alpha=0.25)
strtArray = crossview.plot_array(strt)
cb = plt.colorbar(strtArray, shrink=0.5)
crossview.plot_ibound()

<matplotlib.collections.PatchCollection at 0x1c2b7a39730>

Kh = 1.
Kv = 1.
npf = flopy.mf6.ModflowGwfnpf(gwf, xt3doptions=False,
                                  save_flows=True,
                                  save_specific_discharge=True,
                                  icelltype=1,
                                  k=Kh, k33=Kv)

sto = flopy.mf6.ModflowGwfsto(gwf, sy=0.3, ss=0., iconvert=1)

pd = [(0, 0.7, 0., tModName, 'CONCENTRATION')]
buy = flopy.mf6.ModflowGwfbuy(gwf, denseref=1000., nrhospecies=1,packagedata=pd)

lakeDensity = 1024.5

nlakeconn = 11  # note: number of connections for this lake
# pak_data = [lakeno, strt, nlakeconn, aux, boundname]
pak_data = [(0, 2.25, nlakeconn, 0., 0.)]

connlen = delr / 2.
connwidth = delc
bedleak = 1
con_data = [
    # con_data=(lakeno,iconn,(cellid),claktype,bedleak,belev,telev,connlen,connwidth )
    (0, 0, (0, 0, 0), 'HORIZONTAL', bedleak, 10, 10, connlen, connwidth),
    (0, 1, (1, 0, 1), 'VERTICAL', bedleak, 10, 10, connlen, connwidth),
    (0, 2, (1, 0, 1), 'HORIZONTAL', bedleak, 10, 10, connlen, connwidth),
    (0, 3, (2, 0, 2), 'VERTICAL', bedleak, 10, 10, connlen, connwidth),
    (0, 4, (2, 0, 2), 'HORIZONTAL', bedleak, 10, 10, connlen, connwidth),
    (0, 5, (3, 0, 3), 'VERTICAL', bedleak, 10, 10, connlen, connwidth),
    (0, 6, (2, 0, 4), 'HORIZONTAL', bedleak, 10, 10, connlen, connwidth),
    (0, 7, (2, 0, 4), 'VERTICAL', bedleak, 10, 10, connlen, connwidth),
    (0, 8, (1, 0, 5), 'HORIZONTAL', bedleak, 10, 10, connlen, connwidth),
    (0, 9, (1, 0, 5), 'VERTICAL', bedleak, 10, 10, connlen, connwidth),
    (0, 10, (0, 0, 6), 'HORIZONTAL', bedleak, 10, 10, connlen, connwidth),
]

# period data
p_data = [(0, 'STATUS', 'ACTIVE'), ]

# note: for specifying lake number, use fortran indexing!
fname = '{}.lak.obs.csv'.format(fModName)
lak_obs = {fname: [('lakestage', 'stage', 1),
                   ('lakevolume', 'volume', 1),
                   ('lak1', 'lak', 1, 1),
                   ('lak2', 'lak', 1, 2),
                   ('lak3', 'lak', 1, 3),
                   ('lak4', 'lak', 1, 4),
                   ('lak5', 'lak', 1, 5),
                   ('lak6', 'lak', 1, 6),
                   ('lak7', 'lak', 1, 7),
                   ('lak8', 'lak', 1, 8),
                   ('lak9', 'lak', 1, 9),
                   ('lak10', 'lak', 1, 10),
                   ('lak11', 'lak', 1, 11),
                   ],
           'digits': 10, }

lak = flopy.mf6.modflow.ModflowGwflak(gwf, save_flows=True,
                                      print_input=True, print_flows=True,
                                      print_stage=True,
                                      stage_filerecord='{}.lak.bin'.format(fModName),
                                      budget_filerecord='{}.lak.bud'.format(fModName),
                                      nlakes=len(pak_data),
                                      ntables=0,
                                      packagedata=pak_data, pname='LAK-1',
                                      connectiondata=con_data,
                                      perioddata=p_data,
                                      observations=lak_obs,
                                      auxiliary=['TESTAUXVAR','CONCENTRATION'])

#plot cross section with lake boundary condition
fig, ax = plt.subplots(figsize=(18,6))
crossview = flopy.plot.PlotCrossSection(model=gwf, line={'row':0})
crossview.plot_grid(alpha=0.25)
strtArray = crossview.plot_bc('lak')
cb = plt.colorbar(strtArray, shrink=0.5)
crossview.plot_ibound()

<matplotlib.collections.PatchCollection at 0x1c2b7b3e1c0>

oc = flopy.mf6.ModflowGwfoc(gwf,
                            budget_filerecord='{}.cbc'.format(fModName),
                            head_filerecord='{}.hds'.format(fModName),
                            headprintrecord=[
                                ('COLUMNS', 10, 'WIDTH', 15,
                                 'DIGITS', 6, 'GENERAL')],
                            saverecord=[('HEAD', 'ALL', 'STEPS'),
                                        ('BUDGET', 'ALL', 'STEPS')],
                            printrecord=[('HEAD', 'LAST', 'STEPS'),
                                         ('BUDGET', 'LAST', 'STEPS')])

Transport model (gwt)

#define transport model
gwt = flopy.mf6.ModflowGwt(sim, modelname=tModName)
#define solver
imsgwt = flopy.mf6.ModflowIms(sim, print_option='ALL',
                              outer_dvclose=hclose,
                              outer_maximum=nouter,
                              under_relaxation='NONE',
                              inner_maximum=ninner,
                              inner_dvclose=hclose, rcloserecord=[rclose, 'strict'],
                              linear_acceleration='BICGSTAB',
                              scaling_method='NONE',
                              reordering_method='NONE',
                              relaxation_factor=relax,
                              filename='{}.ims'.format(tModName))
#register solver
sim.register_ims_package(imsgwt, [gwt.name])
#define dis
dis = flopy.mf6.ModflowGwtdis(gwt, nlay=nlay, nrow=nrow, ncol=ncol,
                                      delr=delr, delc=delc,
                                      top=top, botm=botm, idomain=idomain)
# initial conditions
strtConc = 0
ic = flopy.mf6.ModflowGwtic(gwt, strt=strtConc)

# advection
adv = flopy.mf6.ModflowGwtadv(gwt, scheme='UPSTREAM')

# storage
porosity = 0.30
sto = flopy.mf6.ModflowGwtmst(gwt, porosity=porosity)

# sources
sourcerecarray = [(),]
ssm = flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray)

#concentrations from lake.
lak_conc = 35.
lktpackagedata = [(0, lak_conc, 99., 999., 'mylake'), ]
lkt = flopy.mf6.modflow.ModflowGwtlkt(gwt,
                                      boundnames=True,
                                      save_flows=True,
                                      print_input=True,
                                      print_flows=True,
                                      print_concentration=True,
                                      concentration_filerecord=tModName + '.lkt.bin',
                                      budget_filerecord='gwtlak1.bud',
                                      packagedata=lktpackagedata,
                                      pname='LKT-1',
                                      flow_package_name='LAK-1',
                                      flow_package_auxiliary_name='CONCENTRATION',
                                      auxiliary=['aux1', 'aux2'])

# output control
oc = flopy.mf6.ModflowGwtoc(gwt,
                            budget_filerecord='{}.cbc'.format(tModName),
                            concentration_filerecord='{}.ucn'.format(
                                tModName),
                            concentrationprintrecord=[
                                ('COLUMNS', 10, 'WIDTH', 15,
                                 'DIGITS', 6, 'GENERAL')],
                            saverecord=[('CONCENTRATION', 'ALL')],
                            printrecord=[('CONCENTRATION', 'ALL'),
                                         ('BUDGET', 'ALL')])

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

# GWF GWT exchange
gwfgwt = flopy.mf6.ModflowGwfgwt(sim, exgtype='GWF6-GWT6',
                                 exgmnamea=fModName, exgmnameb=tModName,
                                 filename='{}.gwfgwt'.format(name))

Write and run model

sim.write_simulation()
sim.run_simulation()

Model results

This is the code for the model output analysis

import os
import flopy
import numpy as np
import matplotlib.pyplot as plt

name = 'buyLake'
simName = 'buyExample'
fModName = 'buyFlowModel'
tModName = 'buyTrspModel'

ws = os.path.abspath('../Model/')
sim = flopy.mf6.MFSimulation.load(simName,sim_ws=ws)

loading simulation...
  loading simulation name file...
  loading tdis package...
  loading model gwf6...
    loading package dis...
    loading package ic...
    loading package npf...
    loading package sto...
    loading package buy...
    loading package lak...
    loading package oc...
  loading model gwt6...
    loading package dis...
    loading package ic...
    loading package adv...
    loading package mst...
    loading package ssm...
    loading package lkt...
    loading package oc...
    loading package fmi...
  loading exchange package gwf-gwt_exg_0...
  loading ims package buyflowmodel...
  loading ims package buytrspmodel...

fModName, tModName = sim.model_names
print(fModName, tModName)
gwf = sim.get_model(fModName)
gwt = sim.get_model(tModName)

buyflowmodel buytrspmodel

Plot head distributions

headObj = flopy.utils.HeadFile(os.path.join(ws,fModName+'.hds'))
tsList = headObj.get_kstpkper()

def plotHeads(ts):
    tempHead = headObj.get_data(kstpkper=ts)
    fig, ax = plt.subplots(figsize=(18,6))
    crossview = flopy.plot.PlotCrossSection(model=gwf, line={'row':0})
    crossview.plot_grid(alpha=0.25)
    strtArray = crossview.plot_array(tempHead, masked_values=[1e30])
    wtElev = []
    for col in range(tempHead.shape[2]):
        colArray = tempHead[:,0,col]
        wtElev.append(colArray[colArray != 1e30][0])
    wtLine = plt.plot(gwf.modelgrid.xycenters[0],wtElev,c='crimson')
    cb = plt.colorbar(strtArray, shrink=0.5)
    crossview.plot_ibound()

plotHeads(tsList[0])

plotHeads(tsList[49])

Plot concentration distributions

concObj = flopy.utils.HeadFile(os.path.join(ws,tModName+'.ucn'), text='CONCENTRATION')
tsList = concObj.get_kstpkper()

def plotConcs(ts):
    tempConc = concObj.get_data(kstpkper=ts)
    tempHead = headObj.get_data(kstpkper=ts)
    fig, ax = plt.subplots(figsize=(18,6))
    crossview = flopy.plot.PlotCrossSection(model=gwf, line={'row':0})
    crossview.plot_grid(alpha=0.25)
    strtArray = crossview.plot_array(tempConc, masked_values=[1e30])
    wtElev = []
    for col in range(tempHead.shape[2]):
        colArray = tempHead[:,0,col]
        wtElev.append(colArray[colArray != 1e30][0])
    wtLine = plt.plot(gwf.modelgrid.xycenters[0],wtElev,c='crimson')
    cb = plt.colorbar(strtArray, shrink=0.5)
    crossview.plot_ibound()

plotConcs(tsList[0])