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 Pleasant Lake Example; additional examples of working configuration files can be found in the 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 download hydrography first, to ensure that the model area includes all important features. The model should be referenced to a projected coordinate reference system (CRS), ideally with length units of meters and an authority code (such as an EPSG code) 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 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:

 1from mfsetup import MF6model
 2
 3
 4def setup_grid(cfg_file):
 5    """Just set up (a shapefile of) the model grid.
 6    For trying different grid configurations."""
 7    m = MF6model(cfg=cfg_file)
 8    m.setup_grid()
 9    m.modelgrid.write_shapefile('postproc/shps/grid.shp')
10
11if __name__ == '__main__':
12
13    setup_grid('initial_config_poly.yaml')

Download the file: initial_grid_setup.py

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

 1simulation:
 2  sim_name: 'shellmound'
 3  version: 'mf6'
 4  sim_ws: 'model'
 5
 6model:
 7  simulation: 'shellmound'
 8  modelname: 'shellmound'
 9  options:
10    print_input: True
11    save_flows: True
12    newton: True
13  packages: [
14  ]
15
16setup_grid:
17  source_data:
18    features_shapefile:
19      filename: '../mfsetup/tests/data/shellmound/tmr_parent/gis/irregular_boundary.shp'
20  buffer: 0
21  dxy: 1000  # Uniform x, y spacing in meters
22  rotation: 0.
23  crs: 5070  # EPSG code for NAD83 CONUS Albers (meters)
24  snap_to_NHG: True  # option to snap to the USGS National Hydrogeologic Grid

Download the file: 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:

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

Download the file: initial_config_poly.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 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 km2), where NHDPlus HR might have too many lines. NHDPlus version 2 can be obtained from the EPA.

  • NHDPlus High Resolution (HR) 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.

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 (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 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 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 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:

  1. Load all of the tiles to verify that they are correct and cover the whole model active area.

  2. 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).

  1simulation:
  2  sim_name: 'shellmound'
  3  version: 'mf6'
  4  sim_ws: 'model'
  5
  6model:
  7  simulation: 'shellmound'
  8  modelname: 'shellmound'
  9  options:
 10    print_input: True
 11    save_flows: True
 12    newton: True
 13  packages:
 14    - dis
 15    - ic
 16    - np
 17    - oc
 18    - sto
 19    - rch
 20    - sfr
 21    - wel
 22
 23setup_grid:
 24  source_data:
 25    features_shapefile:
 26      filename: '../mfsetup/tests/data/shellmound/tmr_parent/gis/irregular_boundary.shp'
 27  buffer: 0
 28  dxy: 1000  # Uniform x, y spacing in meters
 29  rotation: 0.
 30  crs: 5070  # EPSG code for NAD83 CONUS Albers (meters)
 31  snap_to_NHG: True  # option to snap to the USGS National Hydrogeologic Grid
 32
 33dis:
 34  remake_top: True
 35  options:
 36    length_units: 'meters'
 37  dimensions:
 38    nlay: 1
 39  source_data:
 40    top:
 41      filename: '../mfsetup/tests/data/shellmound/rasters/meras_100m_dem.tif'
 42      elevation_units: 'feet'
 43    botm:
 44      filenames:
 45        0: '../mfsetup/tests/data/shellmound/rasters/mdwy_surf.tif'
 46      elevation_units: 'feet'
 47    idomain:
 48      # polygon shapefile of model active area
 49      filename: '../mfsetup/tests/data/shellmound/tmr_parent/gis/irregular_boundary.shp'
 50
 51tdis:
 52  options:
 53    time_units: 'days'
 54    start_date_time: '2020-01-01'
 55  perioddata:
 56    group 1:
 57      perlen: 1
 58      nper: 1
 59      nstp: 1
 60      steady: True
 61
 62npf:
 63  options:
 64    save_flows: True
 65    rewet: True
 66  griddata:
 67    icelltype: 1
 68    k: 30.
 69    k33: 0.3
 70
 71sto:
 72  options:
 73    save_flows: True
 74  griddata:
 75    iconvert: 1  # convertible layers
 76    sy: 0.2
 77    ss: 1.e-6
 78
 79rch:
 80  options:
 81    print_input: True
 82    print_flows: False
 83    save_flows: True
 84    readasarrays: True
 85  recharge: 0.00025  # 0.00025 m/d ~ 3.5 inches/year
 86
 87sfr:
 88  options:
 89    save_flows: True
 90  source_data:
 91    flowlines:
 92      filename: '../mfsetup/tests/data/shellmound/shps/flowlines.shp'
 93      id_column: 'COMID'  # arguments to sfrmaker.lines.from_shapefile
 94      routing_column: 'tocomid'
 95      width1_column: 'width1'
 96      width2_column: 'width2'
 97      up_elevation_column: 'elevupsmo'
 98      dn_elevation_column: 'elevdnsmo'
 99      name_column: 'GNIS_NAME'
100      width_units: 'feet'  # units of flowline widths
101      elevation_units: 'feet'  # units of flowline elevations
102  sfrmaker_options:
103    one_reach_per_cell: True #  consolidate SFR reaches to one per i, j location
104    to_riv: # convert this line and all downstream lines to the RIV package
105      - 18047206
106
107oc:
108  period_options:
109    0: ['save head last','save budget last']
110
111ims:
112  options:
113    print_option: 'all'
114    complexity: 'complex'
115    csv_output_filerecord: 'solver_out.csv'
116  nonlinear:
117    outer_dvclose: 1.  # m3/d in SFR package
118    outer_maximum: 50
119  linear:
120    inner_maximum: 100
121    inner_dvclose: 0.01
122    rcloserecord: [0.001, 'relative_rclose']

Download the file: 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.

 1import os
 2
 3from mfsetup import MF6model
 4
 5
 6def setup_grid(cfg_file):
 7    """Just set up (a shapefile of) the model grid.
 8    For trying different grid configurations."""
 9    cwd = os.getcwd()
10    m = MF6model(cfg=cfg_file)
11    m.setup_grid()
12    m.modelgrid.write_shapefile('postproc/shps/grid.shp')
13    # Modflow-setup changes the working directory
14    # to the model workspace; change it back
15    os.chdir(cwd)
16
17
18def setup_model(cfg_file):
19    """Set up the whole model."""
20    cwd = os.getcwd()
21    m = MF6model.setup_from_yaml(cfg_file)
22    m.write_input()
23    os.chdir(cwd)
24    return m
25
26
27if __name__ == '__main__':
28
29    #setup_grid('initial_config_poly.yaml')
30    setup_model('initial_config_full.yaml')

Download the file: initial_model_setup.py