08: Modflow-setup demonstration

Modflow-setup is a Python package for rapid, automated construction of MODFLOW models. Often in modeling projects, construction of the basic model structure (discretization, boundary conditions, etc) consumes a lot of time that could be spent on effective history matching, uncertainty quantification and forecast scenario development. While scripting with Flopy can speed construction of the basic model in a way that is robust and repeatable, it still takes time—with each project, new code must be developed or re-worked, dependencies managed, the idiosyncrasies of various interfaces recalled, and inevitable mistakes checked for and debugged.

Modflow-setup aims to speed up construction of the basic model in a way that is robust and repeatable. Grid-independent source data such as shapefiles and rasters are specified in a single configuration file, along with the desired packages and their options. Modflow-setup reads the configuration and the source data, and within a few minutes, produces an external array-based MODFLOW model that is amenable to parameter estimation and uncertainty quantification using PEST (e.g. White et al, 2021). Detailed description of Modflow-setup is provided by Leaf and Fienen (2022) and in online documentation.

Example problem

Here we demonstrate Modflow-setup with a simplified version of the Pleasant Lake model published by Fienen et al (2022), who used a multi-scale modeling approach to evaluate the effects of agricultural groundwater abstraction on the ecology of Pleasant Lake in central Wisconsin, USA. The original report and model files are available at the links below. Figure 1 shows the original published model, which included a coarse (200 meter resolution) “parent” model tightly coupled to a 20 meter resolution local grid refinement (LGR) inset (sub)model (see for example, Mehl and others, 2006; Langevin and others, 2017). The parent model was designed to capture the transient effects of agricultural pumping; the inset model was created to simulate groundwater/lake interactions at a finer scale. To incorporate the effects of distant boundaries, the parent model was in turn loosely coupled to a larger transient regional model via specified head boundaries along the model parameter.

Some advantages of this multi-scale approach over an unstructured grid alternative include:

• The large regional model, which consumes substantial disk space and carries a longer run time, can be managed independently of the parent and inset models. The perimeter of the parent model is sufficiently distant from the area of interest where two-way coupling is not needed.

• Uniform, rectilinear grids are used for the parent and inset models, which greatly simplifies model construction, diagnostics, postprocessing and visualization

• With the parent and inset models tightly coupled at the matrix (inner iteration) level, little or no model solution efficiency is lost vs. an unstructured discretization approach.

Figure 1: The full Pleasant Lake model domain with location map, showing the relationship between the regional, parent and LGR inset models, as well as the irrigation wells considered.

Most details of the Fienen et al (2022) model are included here, but to reduce file sizes and execution time, a smaller parent model domain is used, and the inset model spacing is increased to 40 meters. The example model here is for illustrative purposes only, and is not intended to be adequate for representing pumping impacts or providing a sufficient distance to the perimeter boundaries.

Example model details

• Transient MODFLOW-6 simulation with local grid refinement (LGR)

• “Parent” and LGR “inset” models are dynamically coupled via the Groundwater Flow Exchange Package at the MODFLOW 6 Simulation level

• Coarser “parent” model grid has a uniform 200 meter cell size

• Locally refined “inset” model grid has a uniform 40 meter cell size

• Both parent and inset models are 5 layers

• Parent model is in turn based partially on a larger regional MODFLOW 6 model

  • The RCHa, NPF, STO, and WEL input are rediscretized from the regional model

  • Layer 1 in the regional model is represented by layers 1 and 2 in the parent and inset models

  • The remaining layers are mapped as follows: regional layer 2 to parent/inset layer 3, 3:4 and 4:5

• All models begin with an initial steady-state stress period, followed by monthly timesteps representing calendar year 2012.

• Starting heads for the parent and inset models are resampled from the regional model binary output.

• The Streamflow Routing (SFR) packages in the parent and inset models are constructed from a NHDPlus v2 dataset, and linked together using the Water Mover (MVR) Package

• Head observations are specified from csv files

• LGR inset model extent is based on a buffer distance around a feature of interest (the lake).

• LGR inset model DIS, IC, NPF, STO and RCHa input are rediscretized from the parent model input.

• The Lake package is created from a polygon feature, bathymetry raster, stage-area-volume file and climate data from PRISM.

• Lake package observations are set up automatically (with an output file for each lake)

As we’ll see below, Modflow-setup does all of this automatically, based on input specified in two configuration files (for the parent and inset models) that can be viewed with VSCode or any YAML-aware text editor:

data/pleasant_lgr_parent.yml

data/pleasant_lgr_inset.yml

A final note if you want to try Modflow-setup for your project

The example here is presented in an interactive Jupyter Notebook format primarily for illustrative purposes. Typically, the best way to use Modflow-setup is within a python script (*.py file) that runs from start to finish. To get started with Modflow-setup on your project, see the 10 Minutes to Modflow-setup guide, which includes template configuration files and model build scripts that you can download.

import os
from pathlib import Path
import numpy as np
import pandas as pd
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
from matplotlib import patheffects
import flopy
from mfsetup import MF6model

wd = Path.cwd()

Oftentimes at the start of a modeling project, we want to quickly test different grid resolutions and extents before attempting to build the model. We can do this with Modflow-setup by creating a model instance and then running the setup_grid() method. A model grid instance is created from the setup_grid: block in the configuration file.

%%capture
m = MF6model(cfg='data/pleasant_lgr_parent.yml')
m.setup_grid()

#regional_sim = flopy.mf6.MFSimulation.load('pleasant', sim_ws='../pleasant-lake')
#regional = regional_sim.get_model()
regional = m.parent
regional.modelgrid.write_shapefile('postproc/shps/regional_model_grid.shp')

creating shapely Polygons of grid cells...
finished in 0.07s

writing postproc/shps/regional_model_grid.shp... Done

Since this model has local-grid refinement, it actually consists of two models: a parent built from pleasant_lgr_parent.yml, and an inset built from pleasant_lgr_inset.yml, which is referenced within pleasant_lgr_parent.yml. The two sub-models are connected and solved simulataneously within the same MODFLOW 6 simulation. A model grid is made for each sub-model. The model grids are instances of the MFsetupGrid grid class, a subclass of the Flopy StructuredGrid class with some added functionality.

m.modelgrid

5 layer(s), 25 row(s), 27 column(s)
delr: [200.00...200.00] undefined
delc: [200.00...200.00] undefined
CRS: EPSG:3070
length units: meters
xll: 553200.0; yll: 387800.0; rotation: 0.0
Bounds: (553200.0, 558600.0, 387800.0, 392800.0)

m.inset['plsnt-lgr-inset'].modelgrid

75 row(s), 85 column(s)
delr: [40.00...40.00] undefined
delc: [40.00...40.00] undefined
CRS: EPSG:3070
length units: meters
xll: 554200.0; yll: 388800.0; rotation: 0
Bounds: (554200.0, 557600.0, 388800.0, 391800.0)

Working directory gottcha

Currently, to simplify working with external files in Flopy, Modflow-setup changes the working directory to the model workspace. In the context of a flat script that only builds the model, this is fine, but in a notebook or other workflows, this can potentially cause confusion.

Path.cwd()

PosixPath('/home/runner/work/python-for-hydrology/python-for-hydrology/docs/source/notebooks/part1_flopy/data/pleasant-lgr')

A shapefile of the grid bounding box is written by default on creation of the model grid, to the location specified by output_files: grid_file: in the setup_grid: block (default is <model workspace>/postproc/shps/). A shapefile of the grid cells as polygon features can be written as below:

m.modelgrid.write_shapefile('postproc/shps/plsnt-lgr-parent_grid.shp')
m.inset['plsnt-lgr-inset'].modelgrid.write_shapefile('postproc/shps/plsnt-lgr-inset_grid.shp')

creating shapely Polygons of grid cells...
finished in 0.10s

writing postproc/shps/plsnt-lgr-parent_grid.shp... Done
creating shapely Polygons of grid cells...
finished in 0.10s

writing postproc/shps/plsnt-lgr-inset_grid.shp... Done

Change the working directory back to the notebook location

os.chdir(wd)

Build the whole model

%%capture
m = MF6model.setup_from_yaml('data/pleasant_lgr_parent.yml')

a MF6model instance (subclass of flopy.mf6.ModflowGwf) is returned

m

Pleasant Lake test case version 0.1.post434+ga476e0d
Parent model: /home/runner/work/python-for-hydrology/python-for-hydrology/docs/source/notebooks/part1_flopy/data/pleasant-lake//pleasant
5 layer(s), 25 row(s), 27 column(s)
delr: [200.00...200.00] undefined
delc: [200.00...200.00] undefined
CRS: EPSG:3070
length units: meters
xll: 553200.0; yll: 387800.0; rotation: 0.0
Bounds: (553200.0, 558600.0, 387800.0, 392800.0)
Packages: dis ic npf sto rcha_0 oc chd_0 obs_0 sfr_0 wel_0 obs_1 obs_2
13 period(s):
 per start_datetime end_datetime  perlen  steady  nstp
   0     2012-01-01   2012-01-01     1.0    True     1
   1     2012-01-01   2012-02-01    31.0   False     1
   2     2012-02-01   2012-03-01    29.0   False     1
   ...
  12     2012-12-01   2013-01-01    31.0   False     1

information from the configuration file is stored in an attached cfg dictionary

m.cfg.keys()

dict_keys(['metadata', 'simulation', 'model', 'parent', 'postprocessing', 'setup_grid', 'dis', 'tdis', 'ic', 'npf', 'sto', 'rch', 'sfr', 'high_k_lakes', 'lak', 'mvr', 'chd', 'drn', 'ghb', 'riv', 'wel', 'oc', 'obs', 'ims', 'mfsetup_options', 'filename', 'maw', 'external_files', 'intermediate_data', 'grid'])

the cfg dictionary contains both information from the configuration file, and MODFLOW input (such as external text file arrays) that was developed from the original source data. Internally in Modflow-setup, MODFLOW input in cfg is fed to the various Flopy object constructors.

m.cfg['dis']

defaultdict(dict,
            {'remake_top': True,
             'options': {'length_units': 'meters'},
             'dimensions': {'nlay': 5, 'nrow': 25, 'ncol': 27},
             'griddata': {'top': [{'filename': './external/plsnt-lgr-parent_top.dat'}],
              'botm': [{'filename': './external/plsnt-lgr-parent_botm_000.dat'},
               {'filename': './external/plsnt-lgr-parent_botm_001.dat'},
               {'filename': './external/plsnt-lgr-parent_botm_002.dat'},
               {'filename': './external/plsnt-lgr-parent_botm_003.dat'},
               {'filename': './external/plsnt-lgr-parent_botm_004.dat'}],
              'idomain': [{'filename': './external/plsnt-lgr-parent_idomain_000.dat'},
               {'filename': './external/plsnt-lgr-parent_idomain_001.dat'},
               {'filename': './external/plsnt-lgr-parent_idomain_002.dat'},
               {'filename': './external/plsnt-lgr-parent_idomain_003.dat'},
               {'filename': './external/plsnt-lgr-parent_idomain_004.dat'}]},
             'top_filename_fmt': 'top.dat',
             'botm_filename_fmt': 'botm_{:03d}.dat',
             'idomain_filename_fmt': 'idomain_{:03d}.dat',
             'minimum_layer_thickness': 1,
             'drop_thin_cells': True,
             'source_data': {'top': {'filename': '/home/runner/work/python-for-hydrology/python-for-hydrology/docs/source/notebooks/part1_flopy/data/pleasant-lake/source_data/rasters/dem40m.tif',
               'elevation_units': 'meters'},
              'botm': {'filenames': {1: '/home/runner/work/python-for-hydrology/python-for-hydrology/docs/source/notebooks/part1_flopy/data/pleasant-lake/source_data/rasters/botm0.tif',
                2: '/home/runner/work/python-for-hydrology/python-for-hydrology/docs/source/notebooks/part1_flopy/data/pleasant-lake/source_data/rasters/botm1.tif',
                3: '/home/runner/work/python-for-hydrology/python-for-hydrology/docs/source/notebooks/part1_flopy/data/pleasant-lake/source_data/rasters/botm2.tif',
                4: '/home/runner/work/python-for-hydrology/python-for-hydrology/docs/source/notebooks/part1_flopy/data/pleasant-lake/source_data/rasters/botm3.tif'}}},
             'nlay': 4})

The inset LGR model is attached to the parent within an inset dictionary

m.inset

{'plsnt-lgr-inset': plsnt-lgr-inset model version 0.post434+ga476e0d
 Parent model: ./plsnt-lgr-parent
 5 layer(s), 75 row(s), 85 column(s)
 delr: [40.00...40.00] undefined
 delc: [40.00...40.00] undefined
 CRS: EPSG:3070
 length units: meters
 xll: 554200.0; yll: 388800.0; rotation: 0
 Bounds: (554200.0, 557600.0, 388800.0, 391800.0)
 Packages: dis ic npf sto rcha_0 oc sfr_0 lak_0 obs_0 obs_1
 13 period(s):
  per start_datetime end_datetime  perlen  steady  nstp
    0     2012-01-01   2012-01-01     1.0    True     1
    1     2012-01-01   2012-02-01    31.0   False     1
    2     2012-02-01   2012-03-01    29.0   False     1
    ...
   12     2012-12-01   2013-01-01    31.0   False     1}

Plot the inset and parent model grids with Lake Package connections by layer

inset = m.inset['plsnt-lgr-inset']

# Get the parent model grid extent
# to set plot limits later
l, r, b, t = m.modelgrid.extent

# Make the plot
fig, ax = plt.subplots(figsize=(10, 10))
parent_mv = flopy.plot.PlotMapView(model=m, ax=ax)
inset_mv = flopy.plot.PlotMapView(model=inset, ax=ax)
parent_mv.plot_bc("SFR", plotAll=True)
inset_mv.plot_bc("SFR", plotAll=True)
lcp = parent_mv.plot_grid(lw=0.5, ax=ax)
lci = inset_mv.plot_grid(lw=0.5)

# Get the lake connections
vconn = inset.lak.connectiondata.array[inset.lak.connectiondata.array['claktype'] == 'vertical']
k, i, j = zip(*vconn['cellid'])
lakeconnections = np.zeros((inset.nrow, inset.ncol))
lakeconnections[i, j] = np.array(k)
lakeconnections = np.ma.masked_array(lakeconnections, mask=lakeconnections == 0)
qmi = inset_mv.plot_array(lakeconnections)

# re-limit the plot to the parent model extent
ax.set_ylim(b, t)
ax.set_xlim(l, r)
ax.set_aspect(1)
plt.colorbar(qmi, shrink=0.5, label='Lake connection layer')
ax.set_xlabel('Easting, in WTM meters (epsg: 3070)')
ax.set_ylabel('Northing, in WTM meters (epsg: 3070)')

Text(0, 0.5, 'Northing, in WTM meters (epsg: 3070)')

(just like you would for a Flopy model)

m.write_input()

writing simulation...
  writing simulation name file...
  writing simulation tdis package...
  writing solution package ims_-1...
  writing package pleasant-lgr.mvr...
  writing package pleasant-lgr.gwfgwf...
  writing model plsnt-lgr-parent...
    writing model name file...
    writing package dis...
    writing package ic...
    writing package npf...
    writing package sto...
    writing package rcha_0...
    writing package oc...
    writing package chd_0...
INFORMATION: maxbound in ('gwf6', 'chd', 'dimensions') changed to 272 based on size of stress_period_data
    writing package obs_0...
    writing package sfr_0...
    writing package wel_0...
INFORMATION: maxbound in ('gwf6', 'wel', 'dimensions') changed to 2 based on size of stress_period_data
    writing package obs_1...
    writing package obs_2...
  writing model plsnt-lgr-inset...
    writing model name file...
    writing package dis...
    writing package ic...
    writing package npf...
    writing package sto...
    writing package rcha_0...
    writing package oc...
    writing package sfr_0...
    writing package lak_0...
    writing package obs_0...
    writing package obs_1...
SFRmaker v. 0.11.2

Running Flopy v. 3.7.0 diagnostics...
passed.

Checking for continuity in segment and reach numbering...
passed.

Checking for increasing segment numbers in downstream direction...
passed.

Checking for circular routing...
passed.

Checking reach connections for proximity...
2 segments with non-adjacent reaches found.
At segments:
1 2

2 segments with non-adjacent reaches found.
At segments:
1 2


Checking for model cells with multiple non-zero SFR conductances...
passed.

Checking for streambed tops of less than -10...
passed.

Checking for streambed tops of greater than 15000...
passed.

Checking segment_data for downstream rises in streambed elevation...
Segment elevup and elevdn not specified for nstrm=-17 and isfropt=1
passed.

Checking reach_data for downstream rises in streambed elevation...
7 reaches encountered with strtop < strtop of downstream reach.
Elevation rises:
k i j iseg ireach strtop strtopdn d_strtop reachID diff
4 20 16 2 1 290.66015625 292.6883850097656 2.028228759765625 6 -2.028228759765625
2 24 21 2 10 288.81939697265625 290.66015625 1.84075927734375 15 -1.84075927734375
1 24 22 2 11 288.60980224609375 290.1935119628906 1.583709716796875 16 -1.583709716796875
4 21 16 2 2 290.66015625 292.11285400390625 1.45269775390625 7 -1.45269775390625
3 24 20 2 9 289.3415832519531 290.66015625 1.318572998046875 14 -1.318572998046875
3 24 19 2 8 289.4702453613281 290.66015625 1.189910888671875 13 -1.189910888671875
4 23 19 2 7 290.1935119628906 290.66015625 0.466644287109375 12 -0.466644287109375


Checking reach_data for inconsistencies between streambed elevations and the model grid...
passed.

Checking segment_data for inconsistencies between segment end elevations and the model grid...
Segment elevup and elevdn not specified for nstrm=-17 and isfropt=1
passed.

Checking for streambed slopes of less than 0.0001...
passed.

Checking for streambed slopes of greater than 1.0...
passed.

wrote plsnt-lgr-parent_SFR.chk
wrote plsnt-lgr-parent.sfr.obs
converting reach and segment data to package data...
wrote /home/runner/work/python-for-hydrology/python-for-hydrology/docs/source/notebooks/part1_flopy/data/pleasant-lgr/external/plsnt-lgr-parent_packagedata.dat
wrote /home/runner/work/python-for-hydrology/python-for-hydrology/docs/source/notebooks/part1_flopy/data/pleasant-lgr/plsnt-lgr-parent.sfr

Note: Running the model through Flopy (as below) requires specification of the MODFLOW executable. In Flopy, the executable is specified via the exe_name argument to the simulation constructor for MODFLOW 6, or model constructor for previous MODFLOW versions. Similarly, in Modflow-setup, the exe_name is specified in the simulation: or model: block of the configuration file. This example assumes that a MODFLOW 6 executable with the name “mf6” either resides in the model workspace, or is included in the system path.

m.simulation.run_simulation()

FloPy is using the following executable to run the model: ../../../../../../../../../.local/bin/mf6
                                   MODFLOW 6
                U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
                            VERSION 6.5.0 05/23/2024

   MODFLOW 6 compiled Jun 21 2024 02:57:38 with Intel(R) Fortran Intel(R) 64
   Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0
                             Build 20220726_000000

This software has been approved for release by the U.S. Geological
Survey (USGS). Although the software has been subjected to rigorous
review, the USGS reserves the right to update the software as needed
pursuant to further analysis and review. No warranty, expressed or
implied, is made by the USGS or the U.S. Government as to the
functionality of the software and related material nor shall the
fact of release constitute any such warranty. Furthermore, the
software is released on condition that neither the USGS nor the U.S.
Government shall be held liable for any damages resulting from its
authorized or unauthorized use. Also refer to the USGS Water
Resources Software User Rights Notice for complete use, copyright,
and distribution information.


 MODFLOW runs in SEQUENTIAL mode

 Run start date and time (yyyy/mm/dd hh:mm:ss): 2024/08/07 19:51:54

 Writing simulation list file: mfsim.lst
 Using Simulation name file: mfsim.nam

    Solving:  Stress period:     1    Time step:     1
    Solving:  Stress period:     2    Time step:     1
    Solving:  Stress period:     3    Time step:     1
    Solving:  Stress period:     4    Time step:     1
    Solving:  Stress period:     5    Time step:     1
    Solving:  Stress period:     6    Time step:     1
    Solving:  Stress period:     7    Time step:     1
    Solving:  Stress period:     8    Time step:     1
    Solving:  Stress period:     9    Time step:     1
    Solving:  Stress period:    10    Time step:     1
    Solving:  Stress period:    11    Time step:     1
    Solving:  Stress period:    12    Time step:     1
    Solving:  Stress period:    13    Time step:     1

 Run end date and time (yyyy/mm/dd hh:mm:ss): 2024/08/07 19:51:59
 Elapsed run time:  5.240 Seconds


WARNING REPORT:

  1. OPTIONS BLOCK VARIABLE 'UNIT_CONVERSION' IN FILE 'plsnt-lgr-parent.sfr'
     WAS DEPRECATED IN VERSION 6.4.2. SETTING UNIT_CONVERSION DIRECTLY.
  2. OPTIONS BLOCK VARIABLE 'UNIT_CONVERSION' IN FILE 'plsnt-lgr-inset.sfr'
     WAS DEPRECATED IN VERSION 6.4.2. SETTING UNIT_CONVERSION DIRECTLY.
 Normal termination of simulation.

(True, [])

Path.cwd()

PosixPath('/home/runner/work/python-for-hydrology/python-for-hydrology/docs/source/notebooks/part1_flopy/data/pleasant-lgr')

Read the output and post-process the 3D head results into 2D water table(s)

import flopy.utils.binaryfile as bf
from flopy.utils.postprocessing import get_water_table

lgr_parent_headsobj = bf.HeadFile('plsnt-lgr-parent.hds')
lgr_inset_headsobj = bf.HeadFile('plsnt-lgr-inset.hds')

# read the head results for the last stress period
kper = 12
lgr_parent_hds = lgr_parent_headsobj.get_data(kstpkper=(0, kper))
lgr_inset_hds = lgr_inset_headsobj.get_data(kstpkper=(0, kper))

# Get the water table elevation from the 3D head results
inset_wt = get_water_table(lgr_inset_hds)
parent_wt = get_water_table(lgr_parent_hds)

# put in the lake level (not included in head output)
lake_results = pd.read_csv('lake1.obs.csv')
stage = lake_results['STAGE'][kper]
inset_wt[inset.lakarr[0] == 1] = stage
# add the SFR stage as well
parent_sfr_stage = m.sfr.output.stage().get_data()[0, 0, :]
# get the SFR cell i, j locations
# by unpacking the cellid tuples in the packagedata
psfr_k, psfr_i, psfr_j = zip(*m.sfr.packagedata.array['cellid'])
parent_wt[psfr_i, psfr_j] = parent_sfr_stage
inset_sfr_stage = inset.sfr.output.stage().get_data()[0, 0, :]
# get the SFR cell i, j locations
# by unpacking the cellid tuples in the packagedata
isfr_k, isfr_i, isfr_j = zip(*inset.sfr.packagedata.array['cellid'])
inset_wt[isfr_i, isfr_j] = inset_sfr_stage

# get the cell budget using the .output attribute
# (instead of the binaryfile utility directly)
cbc = m.output.budget()
inset_cbc = inset.output.budget()
lak = inset_cbc.get_data(text='lak', full3D=True)[0].sum(axis=0)
parent_sfr = cbc.get_data(text='sfr', full3D=True)[0]
inset_sfr = inset_cbc.get_data(text='sfr', full3D=True)[0]

Then combine the parent and inset model head results

(into a single grid at the inset model resolution; for a nicer looking plot)

Note: We could skip this step and simply work with the parent and inset models individually in the same plot (as we are doing below with the grids and boundary conditions). However, if we want to create continuous contours of the head solution, we need to first resample it to a single grid (usually at the finer inset model resolution) that spans the whole model domain.

# make the single grid
l, b, r, t = m.modelgrid.bounds
xi = np.arange(l, r, 40)
yi = np.arange(b, t, 40)[::-1]
Xi, Yi = np.meshgrid(xi, yi)

# make a single set of points
# including both parent and inset cell centers
# and water table values
x = m.modelgrid.xcellcenters[~parent_wt.mask]
y = m.modelgrid.ycellcenters[~parent_wt.mask]
x = np.append(x, inset.modelgrid.xcellcenters[~inset_wt.mask])
y = np.append(y, inset.modelgrid.ycellcenters[~inset_wt.mask])
z = parent_wt[~parent_wt.mask].data
z = np.append(z, inset_wt[~inset_wt.mask].data)

# interpolate the results from the points
# onto the single inset resolution grid
results = griddata((x, y), z, (Xi, Yi))

• include the parent and inset model grids

• show the head contours for the combined parent/inset simulation

• show “leakage” results for the Lake and SFR packages (neg. values indicate groundwater discharge to surface water)

Note: Leakage (cell by cell flow) reported by MODFLOW is volumetric, so the apparent parent model values will be higher, simply because they contain more length of stream than the inset model cells.

plt.rcParams['axes.labelsize'] = 10
plt.rcParams['xtick.labelsize'] = 10
plt.rcParams['ytick.labelsize'] = 10

layer = 0
fig, ax = plt.subplots(figsize=(12, 12))
# create Flopy plot objects
parent_mv = flopy.plot.PlotMapView(model=m, ax=ax, layer=layer)
inset_mv = flopy.plot.PlotMapView(model=inset, ax=ax, layer=layer)

vmin, vmax = -100, 100
parent_mv.plot_array(parent_sfr.sum(axis=0), cmap='coolwarm', vmin=vmin, vmax=vmax)
im = inset_mv.plot_array(lak, cmap='coolwarm', vmin=vmin, vmax=vmax)
im = inset_mv.plot_array(inset_sfr.sum(axis=0), cmap='coolwarm', vmin=vmin, vmax=vmax)
cb = fig.colorbar(im, shrink=0.75, label='Leakage, in m$^3$/day')

# contour the combined inset/parent head results
levels = np.arange(290, 315, 2)
ctr = ax.contour(Xi, Yi, results, levels=levels, colors='b', zorder=10)
labels = ax.clabel(ctr, inline=True, fontsize=10, inline_spacing=20)
plt.setp(labels, path_effects=[
    patheffects.withStroke(linewidth=3, foreground="w")])

# plot the grid cell edges
lcp = parent_mv.plot_grid(lw=0.5, ax=ax)
lci = inset_mv.plot_grid(lw=0.5)

ax.set_ylim(b, t)
ax.set_xlim(l, r)
ax.set_aspect(1)
ax.set_ylabel('Northing, in Wisconsin Transverse Mercator (meters)')
ax.set_xlabel('Easting, in Wisconsin Transverse Mercator (meters)')

Text(0.5, 0, 'Easting, in Wisconsin Transverse Mercator (meters)')

References

Fienen, M. N., Haserodt, M. J., Leaf, A. T., and Westenbroek, S. M. (2022). Simulation of regional groundwater flow and groundwater/lake interactions in the central Sands, Wisconsin. U.S. Geological Survey Scientific Investigations Report 2022-5046. doi:10.3133/sir20225046

Langevin, C.D., Hughes, J.D., Banta, E.R., Niswonger, R.G., Panday, S., and Provost, A.M., 2017, Documentation for the MODFLOW 6 groundwater flow model: U.S. Geological Survey Techniques and Methods, book 6, chap. A55, 197 p., https://doi.org/ 10.3133/tm6A55.

Leaf, A.T. and Fienen, M.N. (2022) Modflow-setup: Robust automation of groundwater model construction. Front. Earth Sci. 10:903965. doi: 10.3389/feart.2022.903965

Mehl, S., Hill, M.C., and Leake, S.A., 2006, Comparison of local grid refinement methods for MODFLOW: Groundwater, v. 44, no. 6, p. 792–796, https://doi.org/10.1111/j.1745-6584.2006.00192.x

Westenbroek, S. M., Engott, J. A., Kelson, V. A., and Hunt, R. J. (2018). SWB Version 2.0—a soil-water-balance code for estimating net infiltration and other 1152 water-budget components. U.S. Geological Survey Techniques and Methods, book 6, 118. chap. A59. doi:10.3133/tm6A59

White, J. T., Hemmings, B., Fienen, M. N., and Knowling, M. J. (2021). Towards improved environmental modeling outcomes: Enabling low-cost access to high- 1157 dimensional, geostatistical-based decision-support analyses. Environ. Model. Softw. 139, 105022. doi:10.1016/j.envsoft.2021.105022

Wisconsin Department of Natural Resources (WDNR) (2021). Central Sands Lake study report: Findings and recommendations. Rep. Wis. State Legislature. 1162 doi:10.5281/zenodo.5708791