Specifying boundary conditions with the ‘basic’ MODFLOW stress packages

This page describes configuration file input for the basic MODFLOW stress packages, including the CHD, DRN, GHB, RCH, RIV and WEL packages. The EVT package is not currently supported by Modflow-setup. The supported packages can be broadly placed into two categories. Feature or list-based packages such as CHD, DRN, GHB, RIV and WEL often represent discrete phenomena such as surface water features, pumping wells, or even lines that denote a perimeter boundary. Input to these packages in MODFLOW is tabular, consisting of a table for each stress period, with rows specifying stresses at individual grid cells representing the boundary features. In contrast, continuous or grid-based packages represent a stress field that applies to a large area, such as areal recharge. In past versions of MODFLOW, input to these packages was array-based, with values specified for all model cells, at each stress period. In MODFLOW 6, input to these packages can be array or list-based. The Recharge (RCH) Package is currently the only grid-based stress package supported by Modflow-setup. In keeping with the current structured grid-based paradigm of Modflow-setup, Modflow 6 recharge input is generated for the array-based recharge package (Langevin and others, 2017).

List-based basic stress packages

Input for list-based basic stress packages follows a similar pattern to other packages.

Package blocks are named using the 3 letter MODFLOW abbrieviation for the Package in lower case (e.g. chd:, ghb:, etc.).
Sub-blocks within the package block include:
- options: for specifying Modflow 6 options, exactly as they are described in the input instructions (Langevin and others, 2017).
- source_data: for specifying grid-independent source data to be mapped to the model discretization, in addition to other package input. source_data: in turn can have the following sub-blocks and items:
  
  A shapefile: block for specifying shapefile input that maps the boundary condition features in space. Items in the shapefile block include
  
  filename: path to the shapefile
  
  boundname_col: column in shapefile with feature names to be applied as boundnames in Modflow 6 input
  
  all_touched: argument to rasterio.features.rasterize() that specifies whether all intersected grid cells should be included, or just the grid cells with centers inside the feature.
  
  One or more variable columns: Optionally the shapefile can also supply steady-state variable values by feature in attribute columns named for the variables (e.g. 'head', 'bhead', etc.)
  
  Example:
  
  shapefile: filename: '../../../mfsetup/tests/data/shellmound/shps/waterbodies.shp' id_column: 'COMID' # shapefile attribute field with include_ids: include_ids: [18046162] # features used to form the boundary boundname_column: 'GNIS_NAME'
  
  A csvfile: block for specifying point feature locations or time-varying values for the variables. Items in the shapefile block include
  
  filename: or filenames: path(s) to the csv file(s)
  
  id_column: unique identifier associated with each feature
  
  datetime_column: date-time associated with each stress value
  
  end_datetime_column: date-time associated with the end of stress value (optional; for rates that extend across more than one model stress period. If this is specified, datetime_column: is assumed to indicate the date-time associated with the start of each stress value.
  
  x_col: feature x-coordinate (WEL package only; default 'x')
  
  y_col: feature y-coordinate (WEL package only; default 'y')
  
  length_units: length units associated with the stress value (optional; if omitted no conversion is performed)
  
  time_units: time units associated with the stress value (WEL package only; optional; if omitted no conversion is performed)
  
  volume_units: value-units associated with the stress value (e.g. gallons) in lieu of length-based volume units (e.g., cubic feet) (WEL package only; optional; if omitted volumes are assumed to be in model units of L³ and no conversion is performed)
  
  boundname_col: column in shapefile with feature names to be applied as boundnames in Modflow 6 input
  
  one or more columns for the package variables, specified in the format of <variable>_col, where <variable> is an input variable for the package; for example head_col for the Constant Head Package, or cond_col for the Drain or GHB packages.
  
  period_stats: a sub-block that is used to specify mapping of the input data to the model temporal discretization. Items within period stats are numbered by stress period, with the entry for each item specifying the temporal aggregation. Currently, two options are supported:
  
  aggregation of measurements falling within a stress period. For example, assigning the mean value of all input data points within the stress period. In this case, the aggregration method is simply specified as a string. While mean is typical, any of the standard numpy aggregators can be use (min, max, etc.)
  
  aggregation of measurements from an arbitrary time window. For example, applying a long-term mean to a steady-state stress period, or transient period representing a different time window. In this case three items are specified– the aggregation method, the start date, and end date (e.g. [mean, 2000-01-01, 2017-12-31])
  
  Example:
  
  csvfile: filename: 'shellmound/tables/chd_heads.csv' id_column: 'comid' datetime_column: 'start_datetime' end_datetime_column: 'end_datetime' head_column: 'head' length_units: 'feet' # how heads will be aggregated to the model stress periods period_stats: # apply the mean heads across a specified time period # for the initial steady state # (default of 'mean' uses all times in input data) 0: ['mean', '2000-01-01', '2017-12-31'] 1: mean # apply the mean heads for the time period of s.p. 1 # the last statistic specified (mean) will also be applied to subsequent periods
  
  Additional sub-blocks or items for specifying values for each variable
  
  In general, these sub-blocks are named for the variable (e.g. bhead:).
  
  Scalar values (items) can be specified in model units, and are applied globally to the variable.
  
  Example:
  
  cond: 1.e+3 # in model units
  
  Rasters can be used to specify steady-state values that vary in space; values supplied with a raster are mapped to the model grid using zonal statistics. If the raster contains projection information (GeoTIFFs are preferred in part because of this), reprojection to the model coorindate reference system (CRS) will be performed automatically as needed. Otherwise, the raster is assumed to be in the model projection. Units can optionally be specified and automatically converted; otherwise, the raster values are assumed to be in the model units. Items in the raster block include:
  
  filename: or filenames: path(s) to the raster
  
  length_units: (or elevation_units; optional): length units of the raster values
  
  time_units: (optional): time units of the raster values (cond variable only)
  
  stat: (optional): zonal statistic to use in sampling the raster (defaults are listed for each variable in the Configuration defaults)
  
  Example:
  
  stage: filename: 'shellmound/rasters/meras_100m_dem.tif' elevation_units: 'feet' # zonal statistic to use in sampling elevations in above GeoTIFF to grid stat: 'min' cond: 1.e+3
  
  Not implemented yet: NetCDF input for gridded values that vary in time and space. Due to the lack of standardization in NetCDF coordinate reference information, automatic reprojection is currently not supported for NetCDF files; the data are assumed to be in the model CRS.
- mfsetup_options: Configuration options for Modflow-setup. General options that apply to all basic stress packages include:
  
  external_files: Whether to write the package input as external text arrays or tables (i.e., with open/close statements). By default True except in the case of list-based or tabular files for MODFLOW-NWT models, which are not supported. Adding support for this may require changes to Flopy, which handles external list-based files differently for MODFLOW-2005 style models.
  
  external_filename_fmt: Python string format for external file names. By default, "<package or variable abbreviation>_{:03d}.dat". which results in filenames such as wel_000.dat, wel_001.dat, wel_002.dat… for stress periods 0, 1, and 2, for example.
  
  Other Modflow-setup options specific to individual packages are described below.

Constant Head (CHD) Package

Input consists of specified head values that may vary in time or space.

Required input

parent model head solution –or–

shapefile of features –or–

parent model package (not implemented yet)

at least steady-state head values through one of the methods below

Optional input

raster to specify steady state elevations by cell (for supplied shapefile)

shapefile or csv to specify steady elevations by feature

csv to specify transient elevation by feature (needs to be referenced to features in shapefile)

Examples (also see the Configuration File Gallery)

Setting up a Constant Head package with perimeter fluxes from a parent model (Note: an additional source_data block can be added to represent other features inside of the model perimeter, as below):
chd:
  perimeter_boundary:
    parent_head_file: 'plainfieldlakes/pfl.hds'
Setting up a Constant Head package from features specified in a shapefile, and time-varing heads specified in a csvfile:
chd:
  options:
    save_flows: True
  source_data:
    shapefile:
      filename: '../../../mfsetup/tests/data/shellmound/shps/waterbodies.shp'
      id_column: 'COMID' # shapefile attribute field with include_ids:
      include_ids: [18046162] # features used to form the boundary
      boundname_column: 'GNIS_NAME'
    csvfile:
      filename: 'shellmound/tables/chd_heads.csv'
      id_column: 'comid'
      datetime_column: 'start_datetime'
      end_datetime_column: 'end_datetime'
      head_column: 'head'
      length_units: 'feet'
      # how heads will be aggregated to the model stress periods
      period_stats:
          # apply the mean heads across a specified time period
          # for the initial steady state
          # (default of 'mean' uses all times in input data)
          0: ['mean', '2000-01-01', '2017-12-31']
          1: mean  # apply the mean heads for the time period of s.p. 1
          #  the last statistic specified (mean) will also be applied to subsequent periods

Drain DRN Package

Input consists of elevations and conductances that may vary in time or space.

Required input

shapefile of features –or–

parent model package (not implemented yet)

at least steady-state head and conductance values through one of the methods below

Optional input

global conductance value specified directly

raster to specify steady state elevation by cell (for supplied shapefile)

shapefile or csv to specify steady elevations by feature

csv to specify transient elevation by feature (needs to be referenced to features in shapefile)

Examples (also see the Configuration File Gallery)
drn:
  options:
    save_flows: True
  source_data:
    shapefile:
      filename: '../../../mfsetup/tests/data/shellmound/shps/waterbodies.shp'
      id_column: 'COMID'
      # Include all features in the above shapefile,
      # except those associated with these COMIDs
      exclude_ids: [18046230, 18046226, 18046238, 17953939, 18046140, 18046162]
      boundname_column: 'COMID'
    elev:
      filename: 'shellmound/rasters/meras_100m_dem.tif'
      elevation_units: 'feet'
    cond: 1.e+3  # in model units

General Head Boundary (GHB) Package

Input consists of head elevations and conductances that may vary in time or space.

Required input

shapefile of features –or–
parent model package (not implemented yet)
at least steady-state head and conductance values through one of the methods below

Optional input

global conductance value specified directly
shapefile or csv to specify steady elevations and conductances by feature –or–
rasters to specify steady state elevations or conductances by cell (for supplied shapefile)
csv to specify transient elevations or conductances by feature (needs to be referenced to features in shapefile)

Examples (also see the Configuration File Gallery)

ghb:
  options:
    save_flows: True
  source_data:
    shapefile:
      filename: '../../../mfsetup/tests/data/shellmound/shps/waterbodies.shp'
      id_column: 'COMID'
      # Note: to include all features in a shapefile,
      # simply omit the include_ids: key
      include_ids: [18046230]
      # argument to rasterio.features.rasterize
      # if true, all grid cells touching the feature(s)
      # in the shapefile will be assigned BCs, if False
      # only the cells with centers inside the polygon
      # will be included
      all_touched: True
    # Example of mixed boundary condition input, where
    # a shapefile defines the feature extents
    # a csvfile defines transient head values
    # a raster defines static but spatially varying conductance
    csvfile:
      filename: 'shellmound/tables/chd_heads.csv'
      id_column: 'comid'
      datetime_column: 'start_datetime'
      end_datetime_column: 'end_datetime'
      bhead_column: 'head'
      length_units: 'feet'
      # with no period_stats: block, aggregation defaults
      # to 'mean' for each stress period
    cond:
      # note: a single global value could also be specified here,
      # as for the Drain Package
      filename: shellmound/rasters/k330.tif
      length_units: meters
      time_units: days

River (RIV) package

Input consists of stages, river bottom elevations and conductances,

that may vary in time or space.

Required input

shapefile of features –or–
to_riv: block under sfrmaker_options: with an sfr: block (see configuration gallery)
parent model package (not implemented yet)

Optional input

global conductance value specified directly
default_rbot_thick argument to set a uniform riverbed thickness (rbot = stage - uniform thickness)
shapefile or csv to specify steady heads, conductances and rbots by feature –or–
rasters to specify steady heads, conductances and rbots by cell (for supplied shapefile)
csv to specify transient heads, conductances and rbots by feature (needs to be referenced to features in shapefile)

Examples (also see the Configuration File Gallery)

riv:
  options:
    save_flows: True
  source_data:
    shapefile:
      filename: '../../../mfsetup/tests/data/shellmound/shps/waterbodies.shp'
      id_column: 'COMID'
      include_ids: [17953939]
      # option to include feature names from the shapefile
      # as Modflow boundnames
      boundname_column: 'GNIS_NAME'
    stage:
      filename: 'shellmound/rasters/meras_100m_dem.tif'
      elevation_units: 'feet'
      # zonal statistic to use in sampling elevations in above GeoTIFF to grid
      stat: 'min'
    cond: 1.e+3
  mfsetup_options:
    default_rbot_thickness: 1.

Example of setting up the RIV package using SFRmaker (via the sfr: block):

sfr:
  options:
    save_flows: True
  source_data:
    flowlines:
      filename: '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 source data
      elevation_units: 'feet'  # units of source data
  sfrmaker_options:
    # convert reaches corresponding to these LineString identifiers
    # (in the flowlines id_column), and all downstream reaches
    # to the MODFLOW River package
    to_riv: [18047212]

Well (WEL) Package

Input consists of flux rates that may vary in time or space.

Required input

parent model cell by cell flow solution (not implemented yet) –or–

parent model WEL package

steady-state or transient flux values through one of the methods below

Optional input

temporal discretization (default is to use the average rate(s) for each stress period)

vertical discretization (default is to distribute fluxes vertically by the individual transmissivities of the intersection(s) of the well open interval with the model layers.)

Flux input options with examples (also see the Configuration File Gallery)
Fluxes translated from a parent model WEL package
this input option is very simple. A parent model with a well package is needed, and default_source_data: True must be specified in the parent: block. Then, fluxes from the parent model are simply mapped to the inset model grid, based on the parent model cell centers, and the stress period mappings specified in the parent: block. Well package options can still be specified in a wel: block.
Examples:
wel:
  options:
    print_input: True
    print_flows: True
    save_flows: True
CSV input from one or more files (csvfiles: block)
multiple files can be specified using a list, but column names and units must be consistent

input for column names and units is the same for the general csvfile: block described above

temporal discretization is specified using a period_stats: sub-block

spatial discretization for open intervals spanning multiple layers is specified using a vertical_flux_distribution: sub-block
Examples:
wel:
  options:
    print_input: True
    print_flows: True
    save_flows: True
  source_data:
    csvfiles:
      # pumping input from CSV files
      filenames: ['shellmound/tables/1998-2007_avg_pumping_from_meras21_m3.csv',
                  'shellmound/tables/iwum_m3_6M.csv',
                  'shellmound/tables/sp69_pumping_from_meras21_m3.csv']
      # 'x' and 'y' are the default names for x and y location columns
      volume_units: 'meters'
      time_units: 'days'
      data_column: 'flux_m3'
      datetime_column: 'start_datetime'
      # end datetimes only for input data that needs upsampling;
      # see https://aleaf.github.io/modflow-setup/api/mfsetup.tdis.html#mfsetup.tdis.aggregate_dataframe_to_stress_period
      end_datetime_column: 'end_datetime'
      id_column: 'node'
      period_stats:  # how fluxes will be distributed across model stress periods
        0: none  # no wells simulated in initial period
        1: 'mean'  # mean pumping rate for period 1 and subsequent periods
      vertical_flux_distribution:
        across_layers: False  # False to put fluxes in one layer
        # put wells in layer with thickest or most transmissive intersection with well open interval
        distribute_by: 'transmissivity'  # thickness or transmissivity
        minimum_layer_thickness: 10.  # layers must be at 10 length units thick to have a well;
        # (any dropped wells would be recorded in shellmound_dropped_wells.csv)
Perimeter boundary fluxes from a parent model solution:
wel:
  perimeter_boundary:
    shapefile: 'shellmound/tmr_parent/gis/irregular_boundary.shp'
    parent_cell_budget_file: 'shellmound/tmr_parent/shellmound.cbc'  # needed for the perimeter boundary setup
    parent_binary_grid_file: 'shellmound/tmr_parent/shellmound.dis.grb'
    # parent model head solution
    # for determining boundary fluxes based on saturated thickness
    parent_head_file: 'shellmound/tmr_parent/shellmound.hds'
Similar to the Constant Head Package, a perimeter_boundary block can be mixed with the other input blocks described here to simulate pumping or injection inside of the model perimeter.
wdnr_dataset block
Note

This is a custom option from early versions of Modflow-setup, and is likely to be generalized into a combined shapefile (or CSV site information file) and CSV timeseries input option similar to the other basic stress packages.
site information is specified in a shapefile formatted like csls_sources_wu_pts.shp below

pumping rates are specified by month in a CSV file formatted like master_wu.csv below

temporal discretization is specified with a period_stats: block similar to the csvfiles: option

vertical discretization is specified with a vertical_flux_distribution: block similar to the csvfiles: option
Example:
      water_use: 'plainfieldlakes/source_data/master_wu.csv' # monthly water use rates from WDNR
      water_use_points: 'plainfieldlakes/source_data/wu_points.shp' # point locations of wells in water_use
      period_stats: {0: ['mean', '2012-01-01', '2018-12-31'], # statistic to apply to each stress period 'mean' to average all water use data; <monthname> to average values for a given month across the period (e.g. 'august')
                     1: 'august', # use August pumping rates during test
      }
  output_files:
The vertical_flux_distribution: sub-block

This sub-block specifies how Well Packages fluxes should be distributed vertically.

Items/options include:

across_layers: If True, fluxes for a well will be put in the layer containing the open interval midpoint. If False, fluxes will be distributed to the layers intersecting the well open interval.

distribute_by: 'transmissivity' (default) to distribute fluxes based on the transmissivities of open interval/layer intersections; 'thickness' to distribute fluxes based on intersection thicknesses. Only relevant with across_layers: True.

minimum_layer_thickness: Minimum layer thickness for placing a well (by default 2 model length units). Wells in layers thinner than this will be relocated to the thickess layers at their row, column locations. If no thicker layers exist at the row, column location, the wells are dropped, and reported in <model name>_dropped_wells.csv.

Grid-based basic stress packages

The Recharge (RCH) Package is currently the only grid-based stress package supported by Modflow-setup.

Recharge (RCH) Package

Direct input

As with other grid-based input such as aquifer properties, input to the recharge package can be specified directly as it would in Flopy. This may be useful for setting up a test model quickly. For example, a single scalar value could be entered to apply to all locations across all periods:

rch:
  recharge: 0.001

Or global scalar values could be entered by stress period:

rch:
  recharge:
    0: 0.001
    1: 0.01

In the above example, 0.01 would be also be applied to all subsequent stress periods.

Grid-independent input

Modflow-setup currently supports three methods for entering spatially-referenced recharge input not mapped to the model grid.

Recharge translated from a parent model RCH package

This input option is very simple. A parent model with a recharge package is needed, and default_source_data: True must be specified in the parent: block. Then, fluxes from the parent model are simply mapped to the inset model grid, based on the parent model cell centers, and the stress period mappings specified in the parent: block. Recharge package options can still be specified in a rch: block.
Raster input by stress period
A raster of spatially varying recharge values can be supplied for one or more model stress periods. Similar to the direct input, specified recharge will be applied to subsequent periods were recharge is not specified.

If the raster contains projection information (GeoTIFFs are preferred in part because of this), any reprojection to the model coorindate reference system (CRS) will be performed automatically as needed. Otherwise, the raster is assumed to be in the model projection.

Input items include:

length_units: input recharge length units (optional; if omitted no conversion is performed)

time_units: input recharge time units (optional; if omitted no conversion is performed)

mult: option multiplier value that applies to all stress periods.

resample_method: method for resampling the data from the source grid to model grid. (optional; by default, 'nearest')
Examples:
  source_data:
    rech:
      filenames: # by stress period
        0: 'plainfieldlakes/source_data/net_infiltration__2012-01-01_to_2017-12-31__1066_by_1145__SUM__INCHES_PER_YEAR.tif'
      mult: 0.805
      length_units: 'inches'
      time_units: 'years'
NetCDF input
NetCDF input can be supplied for gridded values that vary in time and space.

Automatic reprojection is supported for Climate Forecast (CF) 1.8-compliant netcdf files (that work with the pyproj.CRS.from_cf() constructor), or files that have a ‘crs_wkt’ or ‘proj4_string’ grid mapping variable (the latter includes many or most Soil Water Balance Code models).

Otherwise, coordinate reference information can be supplied via the crs: item (using any valid input to pyproj.crs.CRS), and the data will be reprojected to the model coordinate reference system.

Input items include:

variable: name of variable in NetCDF file containing the recharge values.

length_units: input recharge length units (optional; if omitted no conversion is performed)

time_units: input recharge time units (optional; if omitted no conversion is performed)

crs: coordinate reference system (CRS) of the netcdf file (optional; only needed if the NetCDF file is in a different CRS than the model and automatic reprojection from the internal grid mapping isn’t working.

resample_method: method for resampling the data from the source grid to model grid. (optional; by default, 'nearest')

period_stats: a sub-block that is used to specify mapping of the input data to the model temporal discretization. Items within period stats are numbered by stress period, with the entry for each item specifying the temporal aggregation. Currently, two options are supported:

aggregation of measurements falling within a stress period. For example, assigning the mean value of all input data points within the stress period. In this case, the aggregration method is simply specified as a string. While mean is typical, any of the standard numpy aggregators can be use (min, max, etc.)

aggregation of measurements from an arbitrary time window. For example, applying a long-term mean to a steady-state stress period, or transient period representing a different time window. In this case three items are specified– the aggregation method, the start date, and end date (e.g. [mean, 2000-01-01, 2017-12-31]; see below for an example)
Examples:
rch:
  options:
    print_input: True
    print_flows: False
    save_flows: True
    readasarrays: True
  source_data:
    # resample recharge from NetCDF file with time-series of grids
    recharge:
      filename: 'shellmound/net_infiltration__2000-01-01_to_2017-12-31__414_by_394.nc'
      variable: 'net_infiltration'
      length_units: 'inches'
      time_units: 'days'
      crs: 5070
      resample_method: 'linear'
      period_stats:
        # for the first two stress periods
        # apply mean recharge rates for dates below
        0: ['mean', '2000-01-01', '2017-12-31']
        1: ['mean', '2000-01-01', '2017-12-31']
        2: 'mean'  # for periods 2 on, use the mean recharge for that period