feat: support for local grid refinement models that only occupy some …

…layers of the enclosing parent model

docs/source/concepts/index.rst

@@ -6,4 +6,5 @@ Modflow-setup concepts and methods
    :maxdepth: 1
 
    Interpolation <interp.rst>
+   Local grid refinement <lgr.rst>
 	Specifying perimeter boundary conditions <perimeter-bcs.rst>

docs/source/concepts/lgr.rst

@@ -0,0 +1,31 @@
+===========================================================
+Local grid refinement
+===========================================================
+
+In MODFLOW 6, two groundwater models can be tightly coupled in the same simulation, which allows for efficient "local grid refinement" (LGR; Mehl and others, 2006) and "semistructured" (Feinstein and others, 2016) configurations that combine multiple structured model layers at different resolutions (Langevin and others, 2017). Locally refined area(s) are conceptualized as separate model(s) that are linked to the surrounding (usually coarser) "parent" model through the GWF Model Exchange Package. Similarly, "semistructured" configurations are represented by multiple linked models for each layer or group of layers with the same resolution.
+
+Modflow-setup supports this functionality via an ``lgr:`` subblock within the ``setup_grid:`` block of the configuration file. The ``lgr:`` subblock consists of one of more subblocks, each keyed by a linked model name and containing configuration input for that linked model. Vertical layer refinement relative to the "parent" model, and the range of parent model layers that will be occupied by the linked "inset" model are also specified.
+
+For example, the following "parent" configuration for the :ref:`Pleasant Lake Example <Pleasant Lake test case>` creates a local refinement model named "``pleasant_lgr_inset``" that spans layers 0 through 3 (the first 4 layers) of the parent model, at the same vertical resolution (1 inset model layer per parent model layer). The horizontal location and resolution of the inset model are specified in the inset model configuration file, in the same way that they are specified for any model.
+
+.. literalinclude:: ../../../examples/pleasant_lgr_parent.yml
+    :language: yaml
+    :start-at: lgr:
+    :end-before: # Structured Discretization Package
+
+Input from the ``lgr:`` subblock and the inset model configuration file(s) is passed to the  :py:class:`Flopy Lgr Utility <flopy.utils.lgrutil.Lgr>`, which helps create input for the GWF Model Exchange Package.
+
+Within the context of a Python session, inset model information is stored in a dictionary under an ``inset`` attribute attached to the parent model. For example, to access a Flopy model object for the above inset model from a parent model named ``model``:
+
+.. code-block:: python
+
+	inset_model = model.inset['pleasant_lgr_inset']
+
+
+**A few notes about LGR functionality in Modflow-setup**
+
+* **Locally refined "inset" models must be aligned with the parent model grid**, which also means that their horizontal resolution must be a factor of the "parent" model resolution. Modflow-setup handles the alignment automatically by "snapping" inset model grids to the lower left corner of the parent cell containing the lower left corner of the inset model grid (the inset model origin in real-world coordinates).
+* Similarly, inset models need to align vertically with the parent model layers. Parent layers can be subdivided using the ``layer_refinement:`` input option.
+* In principal, a **semistructured** configuration consisting of a "parent" model at the coarsest horizontal resolution, and one or more inset models representing [horizontally and potentially vertically] refined layers within the "parent" model domain is possible, but this configuration has not been tested with Modflow-setup. Please contact the Modflow-setup development team by posting a `GitHub issue <https://github.com/DOI-USGS/modflow-setup/issues>`_ if you are interested in this.
+* Similarly, multiple inset models at different horizontal locations, and even inset models within inset models should be achievable, but have not been tested.
+* **Multi-model configurations come with costs.** Each model within a MODFLOW 6 simulation carries its own set of files, including external text array and list input files to packages. As with a single model, when using the automated parameterization functionality in `pyEMU <https://github.com/pypest/pyemu>`_, the number of files is multiplied. At some point, at may be more efficient to work with individual models, and design the grids in such a way that boundary conditions along the model perimeters have a minimal impact on the predictions of interest.

docs/source/conf.py

@@ -236,3 +236,5 @@
     'matplotlib': ('https://matplotlib.org/stable/', None),
     'rasterio': ('https://rasterio.readthedocs.io/en/latest/', None)
 }
+
+tls_verify = False

docs/source/references.rst

@@ -13,6 +13,8 @@ Scripting MODFLOW Model Development Using Python and FloPy: Groundwater, v. 54,
 
 Clark, B.R., Barlow, P.M., Peterson, S.M., Hughes, J.D., Reeves, H.W., and Viger, R.J., 2018, National-scale grid to support regional groundwater availability studies and a national hydrogeologic database: U.S. Geological Survey data release, https://doi.org/10.5066/F7P84B24.
 
+Feinstein, D.T., Fienen, M.N., Reeves, H.W., and Langevin, C.D., 2016, A semi-structured MODFLOW-USG model to evaluate local water sources to wells for decision support: Groundwater, v. 54, no. 4, p. 532–544, accessed June 27, 2017, at https://doi.org/10.1111/gwat.12389.
+
 Fienen, M.N., Haserodt, M.J., Leaf, A.T., Westenbroek, S.M., 2021a, Appendix C: Central sands lakes study technical report: Modeling documentation, Central Sands Lake Study, Wisconsin Department of Natural Resources, 10.5281/zenodo.5708719
 
 Haitjema, H.M. (1995). Analytic Element Modeling of Groundwater Flow. Academic Press, Inc.
@@ -24,6 +26,8 @@ Langevin, C.D., Hughes, J.D., Banta, E.R., Niswonger, R.G., Panday, Sorab, and P
 
 Leaf AT and Fienen MN (2022) Modflow-setup: Robust automation of groundwater model construction. Front. Earth Sci. 10:903965. https://doi.org/10.3389/feart.2022.903965
 
+Mehl, S., Hill, M. C., and Leake, S. A. (2006). Comparison of local grid refinement methods for MODFLOW. Ground Water 44, 792–796. doi:10.1111/j.1745-6584.2006.00192.x
+
 Niswonger, R.G., Panday, S., and Ibaraki, M., 2011, MODFLOW–NWT—A Newton formulation for MODFLOW–2005: U.S. Geological Survey Techniques and Methods, book 6, chap. A37, 44 p. https://doi.org/10.3133/tm6A45
 
 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 water-budget components: U.S. Geological Survey Techniques and Methods, book 6, chap. A59, 118 p., https://doi.org/10.3133/tm6A59.

examples/pleasant_lgr_inset.yml

@@ -44,7 +44,7 @@ dis:
     length_units: 'meters'
   dimensions:
     # nrow and ncol are based on the buffer distance and grid spacing
-    nlay: 5
+    nlay: 4
   source_data:
     top:
       filename: 'data/pleasant/source_data/rasters/dem40m.tif' # DEM file; path relative to setup script
@@ -54,7 +54,6 @@ dis:
         1: 'data/pleasant/source_data/rasters/botm0.tif'  # preprocessed surface for parent model layer 0 bottom
         2: 'data/pleasant/source_data/rasters/botm1.tif'  # preprocessed surface for parent model layer 1 bottom
         3: 'data/pleasant/source_data/rasters/botm2.tif'  # preprocessed surface for parent model layer 2 bottom
-        4: 'data/pleasant/source_data/rasters/botm3.tif'  # preprocessed surface for parent model layer 3 bottom
 
 # Recharge and Well packages are inherited from the parent model
 # (since this is an LGR model

examples/pleasant_lgr_parent.yml

@@ -80,6 +80,8 @@ setup_grid:
     pleasant_lgr_inset:
       filename: 'pleasant_lgr_inset.yml'
       layer_refinement: 1  # number of lgr model layers per parent model layer
+      parent_start_layer: 0 # uppermost parent model layer to include in lgr
+      parent_end_layer: 3 # lowermost parent model layer to include in lgr
 
 # Structured Discretization Package
 dis:

mfsetup/discretization.py

@@ -226,7 +226,7 @@ def get_next_above(seq, value):
 
 
 def fill_cells_vertically(top, botm):
-    """In MODFLOW 6, cells where idomain != 1 are excluded from the solution.
+    """In MODFLOW 6, cells where idomain < 1 are excluded from the solution.
     However, in the botm array, values are needed in overlying cells to
     compute layer thickness (cells with idomain != 1 overlying cells with idomain >= 1 need
     values in botm). Given a 3D numpy array with nan values indicating excluded cells,
@@ -404,14 +404,17 @@ def make_ibound(top, botm, nodata=-9999,
     return idomain
 
 
-def make_lgr_idomain(parent_modelgrid, inset_modelgrid):
+def make_lgr_idomain(parent_modelgrid, inset_modelgrid,
+                     parent_start_layer=0, parent_end_layer=None):
     """Inactivate cells in parent_modelgrid that coincide
     with area of inset_modelgrid."""
     if parent_modelgrid.rotation != inset_modelgrid.rotation:
         raise ValueError('LGR parent and inset models must have same rotation.'
                          f'\nParent rotation: {parent_modelgrid.rotation}'
                          f'\nInset rotation: {inset_modelgrid.rotation}'
                          )
+    if parent_end_layer is None:
+        parent_end_layer = parent_modelgrid.nlay -1
     # upper left corner of inset model in parent model
     # use the cell centers, to avoid edge situation
     # where neighboring parent cell is accidentally selected
@@ -423,7 +426,8 @@ def make_lgr_idomain(parent_modelgrid, inset_modelgrid):
     y1 = inset_modelgrid.ycellcenters[-1, -1]
     pi1, pj1 = parent_modelgrid.intersect(x1, y1, forgive=True)
     idomain = np.ones(parent_modelgrid.shape, dtype=int)
-    idomain[:, pi0:pi1+1, pj0:pj1+1] = 0
+    idomain[parent_start_layer:parent_end_layer,
+            pi0:pi1+1, pj0:pj1+1] = 0
     return idomain
 
 

mfsetup/mf6model.py

@@ -375,11 +375,38 @@ def create_lgr_models(self):
             self.inset[inset_model.name] = inset_model
             #self.inset[inset_model.name]._is_lgr = True
 
-            # create idomain indicating area of parent grid that is LGR
-            lgr_idomain = make_lgr_idomain(self.modelgrid, self.inset[inset_model.name].modelgrid)
-
-            ncpp = int(self.modelgrid.delr[0] / self.inset[inset_model.name].modelgrid.delr[0])
+            # establish inset model layering within parent model
+            parent_start_layer = v.get('parent_start_layer', 0)
+            # parent_end_layer is specified as the last zero-based
+            # parent layer that includes LGR refinement (not as a slice end)
+            parent_end_layer = v.get('parent_end_layer', self.nlay - 1) + 1
             ncppl = v.get('layer_refinement', 1)
+            specified_nlay_lgr = ncppl * (parent_end_layer - parent_start_layer)
+            specified_nlay_dis = inset_cfg['dis']['dimensions']['nlay']
+            if specified_nlay_lgr != specified_nlay_dis:
+                raise ValueError(
+                    f"Parent model layers of {parent_start_layer} to {parent_end_layer} "
+                    f"and layer refinement of {ncppl} implies {specified_nlay_lgr} inset model layers.\n"
+                    f"{specified_nlay_dis} inset model layers specified in DIS package.")
+            # mapping between parent and inset model layers
+            # that is used for copying input from parent model
+            parent_inset_layer_mapping = dict()
+            inset_k = -1
+            for parent_k in range(parent_start_layer, parent_end_layer):
+                for i in range(ncppl):
+                    inset_k += 1
+                    parent_inset_layer_mapping[parent_k] = inset_k
+            self.inset[inset_model.name].cfg['parent']['inset_layer_mapping'] =\
+                parent_inset_layer_mapping
+            # create idomain indicating area of parent grid that is LGR
+            lgr_idomain = make_lgr_idomain(self.modelgrid, self.inset[inset_model.name].modelgrid,
+                                           parent_start_layer, parent_end_layer)
+
+            # inset model horizontal refinement from parent resolution
+            refinement = self.modelgrid.delr[0] / self.inset[inset_model.name].modelgrid.delr[0]
+            if not np.round(refinement, 4).is_integer():
+                raise ValueError(f"LGR inset model spacing must be a factor of the parent model spacing.")
+            ncpp = int(refinement)
             self.lgr[inset_model.name] = Lgr(self.nlay, self.nrow, self.ncol,
                                              self.dis.delr.array, self.dis.delc.array,
                                              self.dis.top.array, self.dis.botm.array,

mfsetup/sourcedata.py

@@ -1425,6 +1425,25 @@ def setup_array(model, package, var, data=None,
                                                 model.cfg['dis']['top_filename_fmt'])
             top = model.load_array(model.cfg['dis']['griddata']['top'][0]['filename'])
 
+        # special case of LGR models
+        # with bottom connections to underlying parent cells
+        if model.version == 'mf6':
+            # (if model is an lgr inset model)
+            if model._is_lgr:
+                # regardless of what is specified for inset model bottom
+                # use top elevations of underlying parent model cells
+                nlay = model.cfg['dis']['dimensions']['nlay']
+                if (nlay < model.parent.modelgrid.nlay):
+                    # use the parent model bottoms
+                    # mapped to the inset model grid by the Flopy Lgr util
+                    lgr = model.parent.lgr[model.name]  # Flopy Lgr inst.
+                    data[nlay-1] = lgr.botm[nlay-1]
+                    # set parent model top to top of
+                    # first active parent model layer below this lgr domain
+                    kp = lgr.get_parent_indices(nlay-1, 0, 0)[0]
+                    model.parent.dis.top = lgr.botmp[kp]
+                    j=2
+
         # fill missing layers if any
         if len(data) < model.nlay:
             all_surfaces = np.zeros((model.nlay + 1, model.nrow, model.ncol), dtype=float) * np.nan

mfsetup/tests/test_lgr.py

@@ -107,8 +107,54 @@ def test_make_lgr_idomain(get_pleasant_lgr_parent_with_grid):
     assert np.all(idomain[:, ~isinset] >= 1)
 
 
-def test_lgr_grid_setup(get_pleasant_lgr_parent_with_grid):
-    m = get_pleasant_lgr_parent_with_grid
+@pytest.mark.parametrize(
+    'lgr_grid_spacing,parent_start_end_layers,inset_nlay', [
+    # parent layers 0 through 3
+    # are specified in pleasant_lgr_parent.yml
+    (40, 'use configuration', 4),
+    (40, (0, 2), 3),
+    # special case to test setup
+    # with no parent start/end layers specified
+    (40, None, 5),
+    pytest.param(35, None, 5, marks=pytest.mark.xfail(
+        reason='inset spacing not a factor of parent spacing')),
+    pytest.param(40, None, 4, marks=pytest.mark.xfail(
+        reason='inset nlay inconsistent with parent layers')),
+    ]
+    )
+def test_lgr_grid_setup(lgr_grid_spacing, parent_start_end_layers,
+                        inset_nlay,
+                        pleasant_lgr_cfg, pleasant_simulation,
+                        project_root_path):
+
+    # apply test parameters to inset/parent configurations
+    inset_cfg_path = project_root_path + '/examples/pleasant_lgr_inset.yml'
+    inset_cfg = load_cfg(inset_cfg_path,
+                         default_file='/mf6_defaults.yml')
+    inset_cfg['setup_grid']['dxy'] = lgr_grid_spacing
+    inset_cfg['dis']['dimensions']['nlay'] = inset_nlay
+
+    cfg = pleasant_lgr_cfg.copy()
+    lgr_cfg = cfg['setup_grid']['lgr']['pleasant_lgr_inset']
+    lgr_cfg['cfg'] = inset_cfg
+    del lgr_cfg['filename']
+    if parent_start_end_layers is None:
+        del lgr_cfg['parent_start_layer']
+        del lgr_cfg['parent_end_layer']
+        k0, k1 = 0, None
+    elif parent_start_end_layers == 'use configuration':
+        k0 = lgr_cfg['parent_start_layer']
+        k1 = lgr_cfg['parent_end_layer']
+    else:
+        k0, k1 = parent_start_end_layers
+        lgr_cfg['parent_start_layer'] = k0
+        lgr_cfg['parent_end_layer'] = k1
+
+    cfg['model']['simulation'] = pleasant_simulation
+    kwargs = get_input_arguments(cfg['model'], mf6.ModflowGwf, exclude='packages')
+    m = MF6model(cfg=cfg, **kwargs)
+    m.setup_grid()
+
     inset_model = m.inset['plsnt_lgr_inset']
     assert isinstance(inset_model, MF6model)
     assert inset_model.parent is m
@@ -120,10 +166,27 @@ def test_lgr_grid_setup(get_pleasant_lgr_parent_with_grid):
         inset_model.modelgrid.write_shapefile(outfolder / 'pleasant_lgr_inset_grid.shp')
 
     # verify that lgr area was removed from parent idomain
-    lgr_idomain = make_lgr_idomain(m.modelgrid, inset_model.modelgrid)
+    lgr_idomain = make_lgr_idomain(m.modelgrid, inset_model.modelgrid,
+                                   k0, k1)
     idomain = m.idomain
     assert idomain[lgr_idomain == 0].sum() == 0
-
+    inset_cells_per_layer = inset_model.modelgrid.shape[1] *\
+        inset_model.modelgrid.shape[2]
+    refinement_factor = int(cfg['setup_grid']['dxy'] / lgr_grid_spacing)
+    nparent_cells_within_lgr_per_layer = inset_cells_per_layer / refinement_factor**2
+    # for each layer where lgr is specified,
+    # there should be at least enough inactive cells to cover the lrg grid area
+    if k1 is None:
+        k1 = m.modelgrid.nlay -1
+    layers_with_lgr = (lgr_idomain == 0).sum(axis=(1, 2)) >=\
+        nparent_cells_within_lgr_per_layer
+    assert all(layers_with_lgr[k0:k1])
+    # layers outside of the specified lgr range should have
+    # less inactive cells than the size of the lgr grid area
+    # (specific to this test case with a large LGR extent relative to the total grid)
+    assert not any(layers_with_lgr[k1+1:])
+    assert not any(layers_with_lgr[:k0])
+    j=2
     # todo: add test that grids are aligned
 
 

mfsetup/wells.py

@@ -83,7 +83,8 @@ def setup_wel_data(model, source_data=None, #for_external_files=True,
 
         # map wells to inset model layers if different from parent
         to_inset_layers = {v:k for k, v in model.parent_layers.items()}
-        spd['k'] = [to_inset_layers[k] for k in spd['k']]
+        spd['k'] = [to_inset_layers.get(k, -9999) for k in spd['k']]
+        spd = spd.loc[spd['k'] >= 0].copy()
 
         df = pd.concat([df, spd], axis=0)