Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

docs: add 10 min to Modflow-setup page #71

Merged
merged 1 commit into from
Aug 22, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
115 changes: 115 additions & 0 deletions docs/source/10min.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,115 @@
10 Minutes to Modflow-setup
============================
This is a short introduction to help get you up and running with Modflow-setup. A complete workflow can be found in the :ref:`Pleasant Lake Example`; additional examples of working configuration files can be found in the :ref:`Configuration File Gallery`.

1) Define the model active area and coordinate reference system
-----------------------------------------------------------------
Depending on the problem, the model area might simply be a box enclosing features of interest and any relevant hydrologic boundaries, or an irregular shape surrounding a watershed or other feature. In either case, it may be helpful to :ref:`download hydrography first <3) Develop flowlines to represent streams>`, to ensure that the model area includes all important features. The model should be referenced to a `projected coordinate reference system (CRS) <https://en.wikipedia.org/wiki/Projected_coordinate_system>`_, ideally with length units of meters and an authority code (such as an `EPSG code <https://en.wikipedia.org/wiki/EPSG_Geodetic_Parameter_Dataset>`_) that unambiguously defines it.

Modflow-setup provides two ways to define a model grid:

* x and y coordinates of the model origin (lower left or upper left corner), grid spacing, number of rows and columns, rotation, and CRS
* As a rectangular area of specified discretization surrounding a polygon shapefile of the model active area (traced by hand or developed by some other means) or a feature of interest buffered by a specified distance.

The active model area is defined subsequently in the DIS package.

.. Note::

Don't forget about the farfield! Usually it is advised to include important competing sinks outside of the immediate area of interest (the nearfield), so that the solution is not over-specified by the perimeter boundary condition, and recognizing that the surface watershed boundary doesn't always coincide exactly with the groundwatershed boundary. See Haitjema (1995) and Anderson and others (2015) for more info.

.. Note::
Need a polygon defining a watershed? In the United States, the `Watershed Boundary Dataset <https://www.usgs.gov/national-hydrography/access-national-hydrography-products>`_ provides watershed deliniations at various scales.


2) Create a setup script and configuration file
------------------------------------------------
Usually creating the desired grid requires some iteration. We can get started on this by making a model setup script and corresponding configuration file.

An initial model setup script for making the model grid:

.. literalinclude:: ../../examples/initial_grid_setup.py
:language: python
:linenos:

Download the file:
:download:`initial_grid_setup.py <../../examples/initial_grid_setup.py>`

An initial configuration file for developing a model grid around a pre-defined active area:

.. literalinclude:: ../../examples/initial_config_poly.yaml
:language: yaml
:linenos:

Download the file:
:download:`initial_config_poly.yaml <../../examples/initial_config_poly.yaml>`

To define a model grid using an origin, grid spacing and dimensions, a ``setup_grid:`` block like this one could be substitued above:

.. literalinclude:: ../../examples/initial_config_box.yaml
:language: yaml
:start-at: setup_grid:

Download the file:
:download:`initial_config_poly.yaml <../../examples/initial_config_box.yaml>`

Now ``initial_setup_script.py`` can be run repeatedly to explore different grids.


3) Develop flowlines to represent streams
------------------------------------------
Next, let's get some data for setting up boundary conditions. For streams, Modflow-setup can accept any linestring shapefile that has a routing column indicating how the lines connect to one another. This can be created by hand, or in the United States, obtained from the National Hydrography Dataset Plus (NHDPlus). There are two types of NHDPlus:

- `NHDPlus version 2 <https://www.epa.gov/waterdata/nhdplus-national-hydrography-dataset-plus>`_ is mapped at the 1:100,000 scale, and is therefore suitable for larger regional models with cell sizes of ~100s of meters to ~1km. NHDPlus version 2 can be the best choice for larger model areas (greater than approx 1,000 km\ :sup:`2`), where NHDPlus HR might have too many lines. NHDPlus version 2 can be obtained from the `EPA <https://www.epa.gov/waterdata/nhdplus-national-hydrography-dataset-plus>`_.
- `NHDPlus High Resolution (HR) <https://www.usgs.gov/national-hydrography/nhdplus-high-resolution>`_ is mapped at the finer 1:24,000 scale, and may therefore work better for smaller problems (discretizations of ~100 meters or less) where better alignment between the mapped lines and stream channel in the DEM is desired, and where the number of linestring features to manage won't be prohibitive. NHDPlus HR can be accessed via the `National Map Downloader <https://apps.nationalmap.gov/downloader/>`_.

Preprocessing NHDPlus HR
^^^^^^^^^^^^^^^^^^^^^^^^^^
Currently, NHDPlus HR data, which comes in a file geodatabase (GDB), must be preprocessed into a shapefile for input to Modflow-setup and `SFRmaker <https://github.com/DOI-USGS/sfrmaker>`_ (which Modflow-setup uses to build the stream network). In many cases, multiple GDBs may need to be combined and undesired line features such as storm sewers culled. The `SFRmaker documentation <https://doi-usgs.github.io/sfrmaker/latest/index.html>`_ has examples for how to read and preprocesses NHDPlus HR.

Preprocessing NHDPlus version 2
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Depending on the application, NHDPlus version 2 may not need to be preprocessed. Reasons to preprocess include:

* the model area is large, and

* read times for one or more NHDPlus drainage basins are slowing the model build
* the DEM being used for the model top is relatively coarse, and sampling a fine DEM during the model build is prohibitive for time or space reasons.

* the stream network is too dense, with too many model cells containing SFR reaches (especially a problem in the eastern US at the 1 km resolution); or there are too many ephemeral streams represented.
* the stream network has divergences where one or more distributary lines are downstream of a confluence.

The `preprocessing module in SFRmaker <https://doi-usgs.github.io/sfrmaker/latest/notebooks/preprocessing_demo.html>`_ can resolve these issues, producing a single set of culled flowlines with width and elevation information and divergences removed. The elevation functionality in the preprocessing module requires a DEM.


4) Get a DEM
-------------
The `National Map Downloader <https://apps.nationalmap.gov/downloader/>`_ has 10 meter DEMs for the United States, with finer resolutions available in many areas. Typically, these come in 1 degree x 1 degree tiles. If many tiles are needed, the uGet Download Manager linked to on the National Map site can automate downloading many tiles. Alternatively, links to the files follow a consistent format, and are therefore amenable to scripted or manual downloads. For example, the tile located between -88 and -87 west and 43 and 44 north is available at:

https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n44w088/USGS_13_n44w088.tif

Making a virtual raster
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Once all of the tiles are downloaded, a virtual raster can be made that allows them to be treated as a single file, without any modifications to the original data. This is required for input to SFRmaker and Modflow-setup. For example, in `QGIS <https://qgis.org/>`_:

a) Load all of the tiles to verify that they are correct and cover the whole model active area.
b) From the ``Raster`` menu, select ``Miscellaneous > Build Virtual Raster``. This will make a virtual raster file with a ``.vrt`` extension that points to the original set of GeoTIFFs, but allows them to be treated as a single continuous raster.

5) Make a minimum working configuration file and model build script
--------------------------------------------------------------------
Now that we have a set of flowlines and a DEM (and perhaps shapefiles for other surface water boundaries), we can fill out the rest of the configuration file to get an initial working model. Later, additional details such as more layers, a well package, observations, or other features can be added in a stepwise approach (Haitjema, 1995).

.. literalinclude:: ../../examples/initial_config_full.yaml
:language: yaml
:linenos:

Download the file:
:download:`initial_config_full.yaml <../../examples/initial_config_full.yaml>`

A setup script for making a minimum working model. Additional functions can be added later to further customize the model outside of the Modflow-setup build step.

.. literalinclude:: ../../examples/initial_model_setup.py
:language: python
:linenos:

Download the file:
:download:`initial_model_setup.py <../../examples/initial_model_setup.py>`
1 change: 1 addition & 0 deletions docs/source/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ Modflow-setup is a Python package for automating the setup of MODFLOW groundwate

Philosophy <philosophy>
Installation <installation>
10 Minutes to Modflow-setup <10min>
Examples <examples>
Configuration File Gallery <config-file-gallery>

Expand Down
3 changes: 3 additions & 0 deletions docs/source/references.rst
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,9 @@ Anderson, M. P., Hunt, R. J., Krohelski, J. T. and Chung, K. (2002),
Using High Hydraulic Conductivity Nodes to Simulate Seepage Lakes. Groundwater, 40: 117-122.
`doi:10.1111/j.1745-6584.2002.tb02496.x <https://doi.org/10.1111/j.1745-6584.2002.tb02496.x>`_

Anderson, M. P., W. W. Woessner, and R. J. Hunt 2015. Applied Groundwater Modeling
620 (Second Edition). San Diego: Academic Press.

Bakker, M., Post, V., Langevin, C. D., Hughes, J. D., White, J. T., Starn, J. J. and Fienen, M. N., 2016,
Scripting MODFLOW Model Development Using Python and FloPy: Groundwater, v. 54, p. 733–739,
`doi:10.1111/gwat.12413. <https://doi.org/10.1111/gwat.12413>`_
Expand Down
24 changes: 24 additions & 0 deletions examples/initial_config_box.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
simulation:
sim_name: 'shellmound'
version: 'mf6'
sim_ws: 'model'

model:
simulation: 'shellmound'
modelname: 'shellmound'
options:
print_input: True
save_flows: True
newton: True
packages: [
]

setup_grid:
xoff: 501405 # lower left x-coordinate
yoff: 1175835 # lower left y-coordinate
nrow: 30
ncol: 35
dxy: 1000
rotation: 0.
epsg: 5070
snap_to_NHG: True
122 changes: 122 additions & 0 deletions examples/initial_config_full.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,122 @@
simulation:
sim_name: 'shellmound'
version: 'mf6'
sim_ws: 'model'

model:
simulation: 'shellmound'
modelname: 'shellmound'
options:
print_input: True
save_flows: True
newton: True
packages:
- dis
- ic
- np
- oc
- sto
- rch
- sfr
- wel

setup_grid:
source_data:
features_shapefile:
filename: '../mfsetup/tests/data/shellmound/tmr_parent/gis/irregular_boundary.shp'
buffer: 0
dxy: 1000 # Uniform x, y spacing in meters
rotation: 0.
crs: 5070 # EPSG code for NAD83 CONUS Albers (meters)
snap_to_NHG: True # option to snap to the USGS National Hydrogeologic Grid

dis:
remake_top: True
options:
length_units: 'meters'
dimensions:
nlay: 1
source_data:
top:
filename: '../mfsetup/tests/data/shellmound/rasters/meras_100m_dem.tif'
elevation_units: 'feet'
botm:
filenames:
0: '../mfsetup/tests/data/shellmound/rasters/mdwy_surf.tif'
elevation_units: 'feet'
idomain:
# polygon shapefile of model active area
filename: '../mfsetup/tests/data/shellmound/tmr_parent/gis/irregular_boundary.shp'

tdis:
options:
time_units: 'days'
start_date_time: '2020-01-01'
perioddata:
group 1:
perlen: 1
nper: 1
nstp: 1
steady: True

npf:
options:
save_flows: True
rewet: True
griddata:
icelltype: 1
k: 30.
k33: 0.3

sto:
options:
save_flows: True
griddata:
iconvert: 1 # convertible layers
sy: 0.2
ss: 1.e-6

rch:
options:
print_input: True
print_flows: False
save_flows: True
readasarrays: True
recharge: 0.00025 # 0.00025 m/d ~ 3.5 inches/year

sfr:
options:
save_flows: True
source_data:
flowlines:
filename: '../mfsetup/tests/data/shellmound/shps/flowlines.shp'
id_column: 'COMID' # arguments to sfrmaker.lines.from_shapefile
routing_column: 'tocomid'
width1_column: 'width1'
width2_column: 'width2'
up_elevation_column: 'elevupsmo'
dn_elevation_column: 'elevdnsmo'
name_column: 'GNIS_NAME'
width_units: 'feet' # units of flowline widths
elevation_units: 'feet' # units of flowline elevations
sfrmaker_options:
one_reach_per_cell: True # consolidate SFR reaches to one per i, j location
to_riv: # convert this line and all downstream lines to the RIV package
- 18047206

oc:
period_options:
0: ['save head last','save budget last']

ims:
options:
print_option: 'all'
complexity: 'complex'
csv_output_filerecord: 'solver_out.csv'
nonlinear:
outer_dvclose: 1. # m3/d in SFR package
outer_maximum: 50
linear:
inner_maximum: 100
inner_dvclose: 0.01
rcloserecord: [0.001, 'relative_rclose']
24 changes: 24 additions & 0 deletions examples/initial_config_poly.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
simulation:
sim_name: 'shellmound'
version: 'mf6'
sim_ws: 'model'

model:
simulation: 'shellmound'
modelname: 'shellmound'
options:
print_input: True
save_flows: True
newton: True
packages: [
]

setup_grid:
source_data:
features_shapefile:
filename: '../mfsetup/tests/data/shellmound/tmr_parent/gis/irregular_boundary.shp'
buffer: 0
dxy: 1000 # Uniform x, y spacing in meters
rotation: 0.
crs: 5070 # EPSG code for NAD83 CONUS Albers (meters)
snap_to_NHG: True # option to snap to the USGS National Hydrogeologic Grid
13 changes: 13 additions & 0 deletions examples/initial_grid_setup.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
from mfsetup import MF6model


def setup_grid(cfg_file):
"""Just set up (a shapefile of) the model grid.
For trying different grid configurations."""
m = MF6model(cfg=cfg_file)
m.setup_grid()
m.modelgrid.write_shapefile('postproc/shps/grid.shp')

if __name__ == '__main__':

setup_grid('initial_config_poly.yaml')
30 changes: 30 additions & 0 deletions examples/initial_model_setup.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
import os

from mfsetup import MF6model


def setup_grid(cfg_file):
"""Just set up (a shapefile of) the model grid.
For trying different grid configurations."""
cwd = os.getcwd()
m = MF6model(cfg=cfg_file)
m.setup_grid()
m.modelgrid.write_shapefile('postproc/shps/grid.shp')
# Modflow-setup changes the working directory
# to the model workspace; change it back
os.chdir(cwd)


def setup_model(cfg_file):
"""Set up the whole model."""
cwd = os.getcwd()
m = MF6model.setup_from_yaml(cfg_file)
m.write_input()
os.chdir(cwd)
return m


if __name__ == '__main__':

#setup_grid('initial_config_poly.yaml')
setup_model('initial_config_full.yaml')