{ "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 automating construction of MODFLOW models from grid-independent source data including shapefiles, rasters, and other MODFLOW models that are geo-located. Input data and model construction options are summarized in a single configuration file in the [YAML](https://en.wikipedia.org/wiki/YAML) format. Source data are read from their native formats and mapped to a regular finite difference grid specified in the configuration file. An external array-based Flopy model instance with the desired packages is created from the sampled source data and configuration settings. Since every modeling problem is different, the Flopy model instance can be further modified in the same script, using Flopy or other Python packages. MODFLOW input can then be written from the flopy model instance.\n", "\n", "## Example problem\n", "Here we use Modflow-setup to make an inset model within the Pleasant Lake example model illustrated in exercise [03: Loading and visualizing groundwater models](solutions/03_Loading_and_visualizing_models-solutions.ipynb). \n", "* The inset model will be the same as the \"parent\" model, except \n", "* the layer top and bottom elevations, head observations and Lake and SFR Package boundary conditions will be re-discretized horizontally on a finer grid. \n", " * The Streamflow Routing (SFR) Packages is constructed from a [NHDPlus v2 dataset](https://nhdplus.com/NHDPlus/NHDPlusV2_data.php)\n", " * The Lake package is created from a polygon feature, bathymetry raster, stage-area-volume file and climate data from [PRISM](https://prismodel.oregonstate.edu/).\n", " * Lake package observations are set up automatically (with an output file for each lake)\n", " * Layer 1 in the parent model is subdivided into layers 1 and 2 in the inset model\n", " * Head observations are (re)specified from csv files\n", "* Array inputs to the RCHa, NPF, STO Packages will be resampled (interpolated) from the parent model to the inset model grid. \n", "* Pumping (Well Package) wells represented in the parent model will be mapped to the nearest cell center in the inset model.\n", "\n", "As we'll see below, Modflow-setup does all of this automatically, based on input specified in a configuration file.\n", "\n", "\n", "**Note:** \n", "The example here is presented in an interactive Jupyter Notebook format 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": "markdown", "metadata": { "vscode": { "languageId": "raw" } }, "source": [ "## Part 1) The configuration file\n", "\n", "The configuration file is the primary vehicle of user input to Modflow-setup. The goal is to have a human and machine-readable document that concisely summarizes the data sources and construction methods/options for the basic model structure prior to parameter estimation. Input is specified in the [YAML](https://en.wikipedia.org/wiki/YAML) format, which is analogous to a Python dictionary but with nesting indicated by indention instead of curly braces. \n", "\n", "* More about YAML [here](https://doi-usgs.github.io/modflow-setup/config-file.html).\n", "* A template configuration file to get started with is available [here](https://doi-usgs.github.io/modflow-setup/10min.html#make-a-minimum-working-configuration-file-and-model-build-script).\n", "\n", "### Contents of `data/pleasant.yml`:\n", "\n", "```yaml\n", "# input for MODFLOW 6 simulation\n", "simulation:\n", " sim_name: 'pleasant-inset'\n", " version: 'mf6'\n", " sim_ws: 'pleasant-inset/'\n", "\n", "# input for MODFLOW 6 model\n", "model:\n", " simulation: 'pleasant-inset'\n", " modelname: 'pleasant-inset'\n", " options:\n", " print_input: True\n", " save_flows: True\n", " newton: True\n", " newton_under_relaxation: True\n", " # packages to build\n", " # (any packages not listed or commented out will not be built,\n", " # event if they have an input block below)\n", " packages: ['dis',\n", " 'ic',\n", " 'npf',\n", " 'oc',\n", " 'sto',\n", " # Note: with no recharge block below and default_source_data=True, \n", " # recharge is regridded from parent model\n", " 'rch', \n", " 'sfr',\n", " 'lak',\n", " 'obs',\n", " 'wel',\n", " 'ims',\n", " 'chd'\n", " ]\n", "\n", "# Regional model to extract boundary conditions, property arrays, and pumping data from\n", "parent:\n", " # arguments to flopy.modflow.Modflow.load for parent model\n", " namefile: 'pleasant.nam'\n", " model_ws: 'pleasant-lake/'\n", " version: 'mf6'\n", " # modflow-setup options\n", " # default_source_data: True, packages and variables that are omitted \n", " # will be pulled from the parent model\n", " default_source_data: True \n", " copy_stress_periods: 'all'\n", " inset_layer_mapping: # mapping between inset and parent model layers\n", " 0: 0 # inset: parent (inset layers 1 and 0 copied from parent layer 0)\n", " 1: 0\n", " 2: 1\n", " 3: 2\n", " 4: 3\n", " start_date_time: '2012-01-01'\n", " length_units: 'meters'\n", " time_units: 'days'\n", " SpatialReference:\n", " xoff: 552400\n", " yoff: 387200\n", " epsg: 3070\n", "\n", "# parameters for setting up the horizontal configuration of the grid\n", "setup_grid:\n", " xoff: 553200 # lower left x-coordinate\n", " yoff: 387800 # lower left y-coordinate\n", " rotation: 0.\n", " dxy: 50 # uniform cell spacing, in CRS units of meters\n", " crs: 3070 # Coordinate reference system (Wisconsin Transverse Mercator)\n", "\n", "# Structured Discretization Package\n", "dis:\n", " options:\n", " length_units: 'meters'\n", " dimensions:\n", " nlay: 5\n", " nrow: 100\n", " ncol: 108\n", " source_data:\n", " top:\n", " # DEM file; path relative to setup script\n", " filename: 'pleasant-lake/source_data/rasters/dem40m.tif'\n", " elevation_units: 'meters'\n", " botm:\n", " filenames: # preprocessed layer surfaces\n", " 1: 'pleasant-lake/source_data/rasters/botm0.tif'\n", " 2: 'pleasant-lake/source_data/rasters/botm1.tif'\n", " 3: 'pleasant-lake/source_data/rasters/botm2.tif'\n", " 4: 'pleasant-lake/source_data/rasters/botm3.tif'\n", "\n", "# Temporal Discretization Package\n", "tdis:\n", " options:\n", " time_units: 'days'\n", " start_date_time: '2012-01-01'\n", " perioddata:\n", " # time discretization info can be specified directly under the perioddata key\n", " # or in groups of stress periods that are discretized in a similar way\n", " group 1: # initial steady-state period (steady specified under sto package)\n", " #perlen: 1 # Specify perlen as an int or list of lengths in model units, \n", " # or perlen=None and 3 of start_date, end_date, nper and/or freq.\"\n", " nper: 1\n", " nstp: 1\n", " tsmult: 1\n", " # \"steady\" can be entered here; \n", " # otherwise the global entry specified in the STO package is used as the default\n", " steady: True \n", " # oc_saverecord: can also be specified by group here; \n", " # otherwise the global entry specified in the oc package is used as the default\n", " group 2: # monthly stress periods\n", " # can be specified by group, \n", " # otherwise start_date_time for the model (under tdis: options) will be used.\n", " start_date_time: '2012-01-01' \n", " # model ends at midnight on this date (2012-12-31 would be the last day simulated)\n", " end_date_time: '2013-01-01'\n", " freq: '1MS' # see arguments to pandas.date_range\n", " nstp: 1\n", " tsmult: 1.5\n", " steady: False\n", "\n", "# Initial Conditions Package\n", "ic:\n", " source_data:\n", " strt:\n", " from_parent: # starting heads sampled from parent model\n", " binaryfile: 'pleasant-lake/pleasant.hds'\n", " stress_period: 0\n", "\n", "# Node Property Flow Package\n", "npf:\n", " options:\n", " save_flows: True\n", " griddata:\n", " icelltype: 1 # variable sat. thickness in all layers\n", " # with parent: default_source_data: True,\n", " # unspecified variables such as \"k\" and \"k33\" are resampled from parent model\n", "\n", "# Storage Package\n", "sto:\n", " options:\n", " save_flows: True\n", " griddata:\n", " iconvert: 1 # convertible layers\n", " # with parent: default_source_data: True,\n", " # unspecified variables such as \"sy\" and \"ss\" are resampled from parent model\n", "\n", "# Well Package\n", "wel:\n", " options:\n", " print_input: True\n", " print_flows: True\n", " save_flows: True\n", " # with parent: default_source_data: True,\n", " # unspecified well fluxes are resampled from parent model\n", "\n", "# Streamflow Routing Package\n", "# SFR input is created using SFRmaker (https://github.com/doi-usgs/sfrmaker)\n", "sfr:\n", " options:\n", " save_flows: True\n", " source_data:\n", " flowlines:\n", " # path to NHDPlus version 2 dataset\n", " nhdplus_paths: ['pleasant-lake/source_data/shps']\n", " # if a DEM is included, streambed top elevations will be sampled\n", " # from the minimum DEM values within a 100 meter buffer around each stream reach\n", " dem:\n", " filename: 'pleasant-lake/source_data/rasters/dem40m.tif'\n", " elevation_units: 'meters'\n", " # SFR observations can be automatically setup from a CSV file\n", " # of x, y locations and fluxes\n", " observations: # see sfrmaker.observations.add_observations for arguments\n", " filename: 'pleasant-lake/source_data/tables/gages.csv'\n", " obstype: 'downstream-flow' # modflow-6 observation type\n", " x_location_column: 'x'\n", " y_location_column: 'y'\n", " obsname_column: 'site_no'\n", " sfrmaker_options:\n", " # the sfrmaker_options: block can include arguments \n", " # to the Lines.to_sfr method in SFRmaker \n", " # or other options such as set_streambed_top_elevations_from_dem\n", " # or to_riv (see SFRmaker and Modflow-setup docs)\n", " set_streambed_top_elevations_from_dem: True\n", "\n", "# Lake Package\n", "lak:\n", " options:\n", " boundnames: True\n", " save_flows: True\n", " surfdep: 0.1 \n", " source_data:\n", " # initial lakebed leakance rates\n", " littoral_leakance: 0.045 # for thin zone around lake perimeter; 1/d\n", " profundal_leakance: 0.025 # for interior of lake basin; 1/d\n", " littoral_zone_buffer_width: 40\n", " lakes_shapefile: # polygon shapefile of lake footprints\n", " filename: 'pleasant-lake/source_data/shps/all_lakes.shp'\n", " id_column: 'HYDROID'\n", " include_ids: [600059060] # pleasant lake\n", " # setup lake precipitation and evaporation from PRISM data\n", " climate:\n", " filenames:\n", " 600059060: 'pleasant-lake/source_data/PRISM_ppt_tmean_stable_4km_189501_201901_43.9850_-89.5522.csv'\n", " format: 'prism'\n", " period_stats:\n", " # for period 0, use average precip and evap for dates below\n", " 0: ['mean', '2012-01-01', '2018-12-31'] # average daily rate for model period for initial steady state\n", " # for subsequent periods,\n", " # average precip and evap to period start/end dates\n", " 1: 'mean' # average daily rate or value for each month\n", " # bathymetry file with lake depths to subtract off model top\n", " bathymetry_raster:\n", " filename: 'pleasant-lake/source_data/rasters/pleasant_bathymetry.tif'\n", " length_units: 'meters'\n", " # bathymetry file with stage/area/volume relationship\n", " stage_area_volume_file:\n", " filename: 'pleasant-lake/source_data/tables/area_stage_vol_Pleasant.csv'\n", " length_units: 'meters'\n", " id_column: 'hydroid'\n", " column_mappings:\n", " volume_m3: 'volume'\n", " external_files: False # option to write connectiondata table to external file\n", "\n", "# Iterative model solution\n", "ims:\n", " options:\n", " complexity: 'complex'\n", " nonlinear:\n", " outer_dvclose: 1.e-1\n", " outer_maximum: 200\n", "\n", "# Observation (OBS) Utility\n", "obs:\n", " source_data:\n", " # Head observations are supplied via csv files \n", " # with x, y locations and observation names\n", " # an observation is generated in each model layer;\n", " # observations at each location can subsequently be processed\n", " # to represent a given screened interval\n", " # for example, using modflow-obs (https://github.com/aleaf/modflow-obs)\n", " filenames: ['pleasant-lake/source_data/tables/nwis_heads_info_file.csv',\n", " 'pleasant-lake/source_data/tables/lake_sites.csv',\n", " 'pleasant-lake/source_data/tables/wdnr_gw_sites.csv',\n", " 'pleasant-lake/source_data/tables/uwsp_heads.csv',\n", " 'pleasant-lake/source_data/tables/wgnhs_head_targets.csv'\n", " ]\n", " column_mappings:\n", " obsname: ['obsprefix', 'obsnme', 'common_name']\n", " drop_observations: ['10019209_lk' # pleasant lake; monitored via Lake Package observations\n", " ]\n", "\n", "# Perimeter boundary condition of heads extracted from parent model solution\n", "chd:\n", " perimeter_boundary:\n", " parent_head_file: 'pleasant-lake/pleasant.hds'\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 2) The model build script\n", "\n", "A template model build script that can be further customized for your project is available [here](https://doi-usgs.github.io/modflow-setup/10min.html#make-a-minimum-working-configuration-file-and-model-build-script).\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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define some functions\n", "\n", "#### A function to 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 Flopy-like model grid instance is created from the ``setup_grid:`` block in the configuration file. \n", "\n", "Wrapping this in a function helps to\n", "1) keep the `__main__` part of the script clean and concise, and \n", "2) also provides a more robust way to control the current working directory (`cwd`). Modflow-setup changes the `cwd` to the specified model workspace (`model_ws:` specified above), to more easily work with external (MODFLOW) input files through the Flopy interface." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def setup_grid(configuration_file):\n", " \"\"\"Just set up (a shapefile of) the model grid.\n", " For trying different grid configurations.\"\"\"\n", " cwd = os.getcwd()\n", " model = MF6model(cfg=configuration_file)\n", " model.setup_grid()\n", " model.modelgrid.write_shapefile('postproc/shps/grid.shp')\n", " os.chdir(cwd)\n", " return model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### A function to build the model using Modflow-setup" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def setup_model(configuration_file):\n", " \"\"\"Set up the whole model.\"\"\"\n", " cwd = os.getcwd()\n", " model = MF6model.setup_from_yaml(configuration_file)\n", " model.write_input()\n", " os.chdir(cwd)\n", " return model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Specify the model configuration file\n", "Read the model workspace location and model name from the configuration file, so that we only have to specify them there." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import yaml\n", "\n", "model_configuration_file = Path('data/pleasant.yml')\n", "# Load the model configuration file into a dictionary\n", "with open(model_configuration_file) as src:\n", " cfg = yaml.load(src, Loader=yaml.Loader)\n", "\n", "sim_ws = model_configuration_file.parent / cfg['simulation']['sim_ws']\n", "model_name = cfg['model']['modelname']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Just make a model grid\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "model = setup_grid('data/pleasant.yml')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(8, 8))\n", "model.modelgrid.plot(ax=ax, zorder=2, ec='b', lw=0.2)\n", "model.parent.modelgrid.plot(ax=ax, lw=0.5)\n", "ax.set_ylim(*model.parent.modelgrid.bounds[1::2])\n", "ax.set_xlim(*model.parent.modelgrid.bounds[::2])\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "parent_model = model.parent\n", "parent_model.modelgrid.write_shapefile(sim_ws / 'postproc/shps/regional_model_grid.shp')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Build the whole model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "model = setup_model('data/pleasant.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": [ "model" ] }, { "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": [ "model.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": [ "model.cfg['dis']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plot the model grid with Lake Package connections by layer" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Get the parent model grid extent\n", "# to set plot limits later\n", "modelgrid = model.modelgrid\n", "l, r, b, t = modelgrid.extent\n", "\n", "# Make the plot\n", "fig, ax = plt.subplots(figsize=(10, 10))\n", "mv = flopy.plot.PlotMapView(model=model, ax=ax)\n", "mv.plot_bc(\"SFR\", plotAll=True)\n", "lcp = mv.plot_grid(lw=0.5, ax=ax)\n", "\n", "# Get the lake connections\n", "vconn = model.lak.connectiondata.array[model.lak.connectiondata.array['claktype'] == 'vertical']\n", "k, i, j = zip(*vconn['cellid'])\n", "lakeconnections = np.zeros((modelgrid.nrow, modelgrid.ncol))\n", "lakeconnections[i, j] = np.array(k) + 1\n", "lakeconnections = np.ma.masked_array(lakeconnections, mask=lakeconnections == 0)\n", "qmi = 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": [ "### Run the model\n", "* we already wrote the input above in the `setup_model()` function (while the working directory was set to the model workspace)\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": [ "model.simulation.run_simulation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the results\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Read the output and post-process the 3D head results into a 2D water table" ] }, { "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", "headfile_object = bf.HeadFile(sim_ws / 'pleasant-inset.hds')\n", "\n", "# read the head results for the last stress period\n", "kper = 12\n", "heads = headfile_object.get_data(kstpkper=(0, kper))\n", "\n", "# Get the water table elevation from the 3D head results\n", "water_table = get_water_table(heads)\n", "\n", "# put in the lake level (not included in head output)\n", "lake_results = pd.read_csv(sim_ws / 'lake1.obs.csv')\n", "stage = lake_results['STAGE'][kper]\n", "cnd = pd.DataFrame(model.lak.connectiondata.array)\n", "k, i, j = zip(*cnd['cellid'])\n", "water_table[i, j] = stage\n", "\n", "# add the SFR stage as well\n", "sfr_stage = model.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", "sfr_k, sfr_i, sfr_j = zip(*model.sfr.packagedata.array['cellid'])\n", "water_table[sfr_i, sfr_j] = sfr_stage\n", "\n", "# get the cell budget using the .output attribute\n", "# (instead of the binaryfile utility directly)\n", "cbc = model.output.budget()\n", "lake_fluxes = cbc.get_data(text='lak', full3D=True)[0].sum(axis=0)\n", "sfr_fluxes = cbc.get_data(text='sfr', full3D=True)[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Make the plot\n", "* show \"leakage\" results for the Lake and SFR packages (neg. values indicate groundwater discharge to surface water)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "levels=np.arange(280, 315, 2)\n", "\n", "fig, ax = plt.subplots(figsize=(8, 8))\n", "pmv = flopy.plot.PlotMapView(model, ax=ax)\n", "ctr = pmv.contour_array(water_table, levels=levels, \n", " linewidths=.5, colors='b')\n", "labels = pmv.ax.clabel(ctr, inline=True, fontsize=8, inline_spacing=1)\n", "vmin, vmax = -100, 100\n", "im = pmv.plot_array(lake_fluxes, cmap='coolwarm', vmin=vmin, vmax=vmax)\n", "im = pmv.plot_array(sfr_fluxes.sum(axis=0), cmap='coolwarm', vmin=vmin, vmax=vmax)\n", "cb = fig.colorbar(im, shrink=0.7, label='Leakage, in m$^3$/day')\n", "ax.set_ylabel(\"Northing, WTM meters\")\n", "ax.set_xlabel(\"Easting, WTM meters\")\n", "ax.set_aspect(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What if we want to further customize the model?\n", "\n", "**Add a high-capacity pumping well at a specified location**, using the Multi-Aquifer Well (MAW) Package, which isn't yet supported by Modflow-setup" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "well_x, well_y = 556000, 389500\n", "# get the row, column location from the coordinates\n", "well_i, well_j = model.modelgrid.intersect(well_x, well_y)\n", "well_i, well_j" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set up the MAW Package" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#wellno radius bottom strt condeqn ngwfnodes boundname\n", "packagedata = [[0, 0.25, 274, 296., 'mean', 4, 'new_well']]\n", "#wellno icon k i j scrn_top scrn_botm hk_skin radius_skin\n", "connectiondata = [\n", " [0, 0, 0, well_i, well_j, 299, 274, 100.0, 0.7],\n", " [0, 1, 1, well_i, well_j, 299, 274, 100.0, 0.7],\n", " [0, 2, 2, well_i, well_j, 299, 274, 100.0, 0.7],\n", " [0, 3, 3, well_i, well_j, 299, 274, 100.0, 0.7]\n", " ]\n", "perioddata_input = [[0, \"rate\", -2000*1440/264.17]]\n", "\n", "maw = flopy.mf6.ModflowGwfmaw(\n", " model,\n", " boundnames=True,\n", " save_flows=True,\n", " packagedata=packagedata,\n", " connectiondata=connectiondata,\n", " perioddata=perioddata_input,\n", " pname='maw'\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Write the MAW Package input file;\n", "\n", "Refresh the Name file to include the MAW Package" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model.maw.write()\n", "# refresh the name file with the added package\n", "model.name_file.write()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model.simulation.run_simulation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the simulated drawdown" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "heads2 = headfile_object.get_data(kstpkper=(0, kper))\n", "\n", "\n", "water_table2 = get_water_table(heads2)\n", "lake_results = pd.read_csv(sim_ws / 'lake1.obs.csv')\n", "stage = lake_results['STAGE'][kper]\n", "cnd = pd.DataFrame(model.lak.connectiondata.array)\n", "k, i, j = zip(*cnd['cellid'])\n", "water_table2[i, j] = stage\n", "sfr_stage = model.sfr.output.stage().get_data()[0, 0, :]\n", "sfr_k, sfr_i, sfr_j = zip(*model.sfr.packagedata.array['cellid'])\n", "water_table2[sfr_i, sfr_j] = sfr_stage\n", "drawdown = water_table - water_table2\n", "\n", "# get the cell budget using the .output attribute\n", "# (instead of the binaryfile utility directly)\n", "cbc = model.output.budget()\n", "lake_fluxes = cbc.get_data(text='lak', full3D=True)[0].sum(axis=0)\n", "sfr_fluxes = cbc.get_data(text='sfr', full3D=True)[0]\n", "\n", "levels=np.arange(1, 10 , 1)\n", "\n", "fig, ax = plt.subplots(figsize=(8, 8))\n", "pmv = flopy.plot.PlotMapView(model, ax=ax)\n", "ctr = pmv.contour_array(drawdown, levels=levels, \n", " linewidths=.5, colors='b')\n", "labels = pmv.ax.clabel(ctr, inline=True, fontsize=8, inline_spacing=1)\n", "vmin, vmax = -100, 100\n", "im = pmv.plot_array(lake_fluxes, cmap='coolwarm', vmin=vmin, vmax=vmax)\n", "im = pmv.plot_array(sfr_fluxes.sum(axis=0), cmap='coolwarm', vmin=vmin, vmax=vmax)\n", "cb = fig.colorbar(im, shrink=0.7, label='Leakage, in m$^3$/day')\n", "ax.set_ylabel(\"Northing, WTM meters\")\n", "ax.set_xlabel(\"Easting, WTM meters\")\n", "ax.set_aspect(1)" ] } ], "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.12.11" } }, "nbformat": 4, "nbformat_minor": 4 }