{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 08: Modflow-setup demonstration\n", "\n", "[Modflow-setup](https://github.com/DOI-USGS/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. \n", "\n", "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](https://doi.org/10.1016/j.envsoft.2021.105022)). Detailed description of Modflow-setup is provided by [Leaf and Fienen (2022)](https://www.frontiersin.org/articles/10.3389/feart.2022.903965) and in [online documentation](https://doi-usgs.github.io/modflow-setup/latest/).\n", "\n", "#### Example problem\n", "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. \n", "\n", "##### Some advantages of this multi-scale approach over an unstructured grid alternative include:\n", "* 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.\n", "* Uniform, rectilinear grids are used for the parent and inset models, which greatly simplifies model construction, diagnostics, postprocessing and visualization\n", "* 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.\n", "\n", "\n", "\n", "**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.\n", "\n", "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.\n", "\n", "##### Example model details \n", "\n", "* Transient MODFLOW-6 simulation with local grid refinement (LGR)\n", "* \"Parent\" and LGR \"inset\" models are dynamically coupled via the Groundwater Flow Exchange Package at the MODFLOW 6 Simulation level\n", "* Coarser \"parent\" model grid has a uniform 200 meter cell size \n", "* Locally refined \"inset\" model grid has a uniform 40 meter cell size\n", "* Both parent and inset models are 5 layers\n", "* Parent model is in turn based partially on a larger regional MODFLOW 6 model\n", " * The RCHa, NPF, STO, and WEL input are rediscretized from the regional model\n", " * Layer 1 in the regional model is represented by layers 1 and 2 in the parent and inset models\n", " * The remaining layers are mapped as follows: regional layer 2 to parent/inset layer 3, 3:4 and 4:5\n", "* All models begin with an initial steady-state stress period, followed by monthly timesteps representing calendar year 2012.\n", "* Starting heads for the parent and inset models are resampled from the regional model binary output.\n", "* The Streamflow Routing (SFR) packages in the parent and inset models are constructed from a [NHDPlus v2 dataset](https://nhdplus.com/NHDPlus/NHDPlusV2_data.php), and linked together using the Water Mover (MVR) Package\n", "* Head observations are specified from csv files\n", "* LGR inset model extent is based on a buffer distance around a feature of interest (the lake).\n", "* LGR inset model DIS, IC, NPF, STO and RCHa input are rediscretized from the parent model input.\n", "* The Lake package is created from a polygon feature, bathymetry raster, stage-area-volume file and climate data from [PRISM](https://prism.oregonstate.edu/).\n", "* Lake package observations are set up automatically (with an output file for each lake)\n", "\n", "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:\n", " \n", "``data/pleasant_lgr_parent.yml`` \n", "``data/pleasant_lgr_inset.yml`` \n", "\n", "\n", "#### A final note if you want to try Modflow-setup for your project\n", "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](https://doi-usgs.github.io/modflow-setup/latest/10min.html) guide, which includes template configuration files and model build scripts that you can download.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "from pathlib import Path\n", "import numpy as np\n", "import pandas as pd\n", "from scipy.interpolate import griddata\n", "import matplotlib.pyplot as plt\n", "from matplotlib import patheffects\n", "import flopy\n", "from mfsetup import MF6model\n", "\n", "wd = Path.cwd()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Just make a model grid\n", "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. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "m = MF6model(cfg='data/pleasant_lgr_parent.yml')\n", "m.setup_grid()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#regional_sim = flopy.mf6.MFSimulation.load('pleasant', sim_ws='../pleasant-lake')\n", "#regional = regional_sim.get_model()\n", "regional = m.parent\n", "regional.modelgrid.write_shapefile('postproc/shps/regional_model_grid.shp')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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 simultaneously 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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m.modelgrid" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m.inset['plsnt-lgr-inset'].modelgrid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Working directory gottcha\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Path.cwd()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Write shapefiles of the inset and parent modelgrids\n", "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 ``/postproc/shps/``). A shapefile of the grid cells as polygon features can be written as below:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m.modelgrid.write_shapefile('postproc/shps/plsnt-lgr-parent_grid.shp')\n", "m.inset['plsnt-lgr-inset'].modelgrid.write_shapefile('postproc/shps/plsnt-lgr-inset_grid.shp')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Change the working directory back to the notebook location" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "os.chdir(wd)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Build the whole model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "m = MF6model.setup_from_yaml('data/pleasant_lgr_parent.yml')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "a ``MF6model`` instance (subclass of ``flopy.mf6.ModflowGwf``) is returned" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "information from the configuration file is stored in an attached ``cfg`` dictionary" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m.cfg.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m.cfg['dis']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The inset LGR model is attached to the parent within an ``inset`` dictionary" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m.inset" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plot the inset and parent model grids with Lake Package connections by layer" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inset = m.inset['plsnt-lgr-inset']\n", "\n", "# Get the parent model grid extent\n", "# to set plot limits later\n", "l, r, b, t = m.modelgrid.extent\n", "\n", "# Make the plot\n", "fig, ax = plt.subplots(figsize=(10, 10))\n", "parent_mv = flopy.plot.PlotMapView(model=m, ax=ax)\n", "inset_mv = flopy.plot.PlotMapView(model=inset, ax=ax)\n", "parent_mv.plot_bc(\"SFR\", plotAll=True)\n", "inset_mv.plot_bc(\"SFR\", plotAll=True)\n", "lcp = parent_mv.plot_grid(lw=0.5, ax=ax)\n", "lci = inset_mv.plot_grid(lw=0.5)\n", "\n", "# Get the lake connections\n", "vconn = inset.lak.connectiondata.array[inset.lak.connectiondata.array['claktype'] == 'vertical']\n", "k, i, j = zip(*vconn['cellid'])\n", "lakeconnections = np.zeros((inset.nrow, inset.ncol))\n", "lakeconnections[i, j] = np.array(k)\n", "lakeconnections = np.ma.masked_array(lakeconnections, mask=lakeconnections == 0)\n", "qmi = inset_mv.plot_array(lakeconnections)\n", "\n", "# re-limit the plot to the parent model extent\n", "ax.set_ylim(b, t)\n", "ax.set_xlim(l, r)\n", "ax.set_aspect(1)\n", "plt.colorbar(qmi, shrink=0.5, label='Lake connection layer')\n", "ax.set_xlabel('Easting, in WTM meters (epsg: 3070)')\n", "ax.set_ylabel('Northing, in WTM meters (epsg: 3070)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### write the MODFLOW input files\n", "(just like you would for a Flopy model)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m.write_input()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run the model\n", "\n", "**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](https://doi-usgs.github.io/modflow-setup/latest/config-file-gallery.html#pleasant-lake-test-case). 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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m.simulation.run_simulation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the results\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Path.cwd()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Read the output and post-process the 3D head results into 2D water table(s)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import flopy.utils.binaryfile as bf\n", "from flopy.utils.postprocessing import get_water_table\n", "\n", "lgr_parent_headsobj = bf.HeadFile('plsnt-lgr-parent.hds')\n", "lgr_inset_headsobj = bf.HeadFile('plsnt-lgr-inset.hds')\n", "\n", "# read the head results for the last stress period\n", "kper = 12\n", "lgr_parent_hds = lgr_parent_headsobj.get_data(kstpkper=(0, kper))\n", "lgr_inset_hds = lgr_inset_headsobj.get_data(kstpkper=(0, kper))\n", "\n", "# Get the water table elevation from the 3D head results\n", "inset_wt = get_water_table(lgr_inset_hds)\n", "parent_wt = get_water_table(lgr_parent_hds)\n", "\n", "# put in the lake level (not included in head output)\n", "lake_results = pd.read_csv('lake1.obs.csv')\n", "stage = lake_results['STAGE'][kper]\n", "inset_wt[inset.lakarr[0] == 1] = stage\n", "# add the SFR stage as well\n", "parent_sfr_stage = m.sfr.output.stage().get_data()[0, 0, :]\n", "# get the SFR cell i, j locations\n", "# by unpacking the cellid tuples in the packagedata\n", "psfr_k, psfr_i, psfr_j = zip(*m.sfr.packagedata.array['cellid'])\n", "parent_wt[psfr_i, psfr_j] = parent_sfr_stage\n", "inset_sfr_stage = inset.sfr.output.stage().get_data()[0, 0, :]\n", "# get the SFR cell i, j locations\n", "# by unpacking the cellid tuples in the packagedata\n", "isfr_k, isfr_i, isfr_j = zip(*inset.sfr.packagedata.array['cellid'])\n", "inset_wt[isfr_i, isfr_j] = inset_sfr_stage\n", "\n", "# get the cell budget using the .output attribute\n", "# (instead of the binaryfile utility directly)\n", "cbc = m.output.budget()\n", "inset_cbc = inset.output.budget()\n", "lak = inset_cbc.get_data(text='lak', full3D=True)[0].sum(axis=0)\n", "parent_sfr = cbc.get_data(text='sfr', full3D=True)[0]\n", "inset_sfr = inset_cbc.get_data(text='sfr', full3D=True)[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Then combine the parent and inset model head results\n", "(into a single grid at the inset model resolution; for a nicer looking plot)\n", "\n", "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. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# make the single grid\n", "l, b, r, t = m.modelgrid.bounds\n", "xi = np.arange(l, r, 40)\n", "yi = np.arange(b, t, 40)[::-1]\n", "Xi, Yi = np.meshgrid(xi, yi)\n", "\n", "# make a single set of points\n", "# including both parent and inset cell centers\n", "# and water table values\n", "x = m.modelgrid.xcellcenters[~parent_wt.mask]\n", "y = m.modelgrid.ycellcenters[~parent_wt.mask]\n", "x = np.append(x, inset.modelgrid.xcellcenters[~inset_wt.mask])\n", "y = np.append(y, inset.modelgrid.ycellcenters[~inset_wt.mask])\n", "z = parent_wt[~parent_wt.mask].data\n", "z = np.append(z, inset_wt[~inset_wt.mask].data)\n", "\n", "# interpolate the results from the points\n", "# onto the single inset resolution grid\n", "results = griddata((x, y), z, (Xi, Yi))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Make the plot\n", "* include the parent and inset model grids\n", "* show the head contours for the combined parent/inset simulation\n", "* show \"leakage\" results for the Lake and SFR packages (neg. values indicate groundwater discharge to surface water)\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['axes.labelsize'] = 10\n", "plt.rcParams['xtick.labelsize'] = 10\n", "plt.rcParams['ytick.labelsize'] = 10\n", "\n", "layer = 0\n", "fig, ax = plt.subplots(figsize=(12, 12))\n", "# create Flopy plot objects\n", "parent_mv = flopy.plot.PlotMapView(model=m, ax=ax, layer=layer)\n", "inset_mv = flopy.plot.PlotMapView(model=inset, ax=ax, layer=layer)\n", "\n", "vmin, vmax = -100, 100\n", "parent_mv.plot_array(parent_sfr.sum(axis=0), cmap='coolwarm', vmin=vmin, vmax=vmax)\n", "im = inset_mv.plot_array(lak, cmap='coolwarm', vmin=vmin, vmax=vmax)\n", "im = inset_mv.plot_array(inset_sfr.sum(axis=0), cmap='coolwarm', vmin=vmin, vmax=vmax)\n", "cb = fig.colorbar(im, shrink=0.75, label='Leakage, in m$^3$/day')\n", "\n", "# contour the combined inset/parent head results\n", "levels = np.arange(290, 315, 2)\n", "ctr = ax.contour(Xi, Yi, results, levels=levels, colors='b', zorder=10)\n", "labels = ax.clabel(ctr, inline=True, fontsize=10, inline_spacing=20)\n", "plt.setp(labels, path_effects=[\n", " patheffects.withStroke(linewidth=3, foreground=\"w\")])\n", "\n", "# plot the grid cell edges\n", "lcp = parent_mv.plot_grid(lw=0.5, ax=ax)\n", "lci = inset_mv.plot_grid(lw=0.5)\n", "\n", "ax.set_ylim(b, t)\n", "ax.set_xlim(l, r)\n", "ax.set_aspect(1)\n", "ax.set_ylabel('Northing, in Wisconsin Transverse Mercator (meters)')\n", "ax.set_xlabel('Easting, in Wisconsin Transverse Mercator (meters)')\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### References\n", "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\n", "\n", "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.\n", "\n", "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\n", "\n", "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\n", "\n", "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,\n", "118. chap. A59. doi:10.3133/tm6A59\n", "\n", "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.\n", "139, 105022. doi:10.1016/j.envsoft.2021.105022\n", "\n", "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" ] } ], "metadata": { "kernelspec": { "display_name": "pyclass", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.7" } }, "nbformat": 4, "nbformat_minor": 4 }