.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/mf6/hondsrug.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_mf6_hondsrug.py: Regional model ============== This example shows a simplified script for building a groundwater model in the northeast of the Netherlands. A primary feature of this area is an ice-pushed ridge called the Hondsrug. This examples demonstrates modifying external data for use in a MODFLOW6 model. In overview, the model features: * Thirteen layers: seven aquifers and six aquitards; * A dense ditch network in the east; * Pipe drainage for agriculture; * Precipitation and evapotranspiration. .. GENERATED FROM PYTHON SOURCE LINES 18-20 .. code-block:: default .. GENERATED FROM PYTHON SOURCE LINES 22-23 We'll start with the usual imports, and an import from scipy. .. GENERATED FROM PYTHON SOURCE LINES 23-32 .. code-block:: default import matplotlib.pyplot as plt import numpy as np import pandas as pd import scipy.ndimage import xarray as xr import imod .. GENERATED FROM PYTHON SOURCE LINES 33-37 Before starting to create the input data, we will create the groundwater model variable (gwf_model) using `imod.mf6.GroundwaterFlowModel `_. The data from all the model packages will be added to this variable. .. GENERATED FROM PYTHON SOURCE LINES 37-40 .. code-block:: default gwf_model = imod.mf6.GroundwaterFlowModel() .. GENERATED FROM PYTHON SOURCE LINES 41-48 This package allows specifying a regular MODFLOW grid. This grid is assumed to be rectangular horizontally, but can be distorted vertically. Load data --------- We'll load the data from the examples that come with this package. .. GENERATED FROM PYTHON SOURCE LINES 48-62 .. code-block:: default layermodel = imod.data.hondsrug_layermodel() # Make sure that the idomain is provided as integers idomain = layermodel["idomain"].astype(int) # We only need to provide the data for the top as a 2D array. Modflow 6 will # compare the top against the uppermost active bottom cell. top = layermodel["top"].max(dim="layer") bot = layermodel["bottom"] top.plot.imshow() .. image-sg:: /examples/mf6/images/sphx_glr_hondsrug_001.png :alt: dx = 25.0, dy = -25.0 :srcset: /examples/mf6/images/sphx_glr_hondsrug_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 63-74 Adding information to the DIS package ------------------------------------- The following step is to add the previously created discretization data to the gwf_model variable. This is done using the function `imod.mf6.StructuredDiscretization `_. The data to include is the top of the model domain, the bottom of the layers, and the idomain. All this information comes from the previously imported tifs (now converted to `xarray.DataArray `_. .. GENERATED FROM PYTHON SOURCE LINES 74-79 .. code-block:: default gwf_model["dis"] = imod.mf6.StructuredDiscretization( top=top, bottom=bot, idomain=idomain ) .. GENERATED FROM PYTHON SOURCE LINES 80-90 Node property flow package - NPF ================================= This package contains the information related to the aquifer properties used to calculate hydraulic conductance. This package replaces the Layer Property Flow (LPF), Block-Centered Flow (BCF), and Upstream Weighting (UPW) packages from previous MODFLOW versions. Hydraulic conductivity ---------------------- .. GENERATED FROM PYTHON SOURCE LINES 90-92 .. code-block:: default k = layermodel["k"] .. GENERATED FROM PYTHON SOURCE LINES 93-119 icelltype ---------- The cell type to be used in the model (confined or convertible) can be defined under ICELLTYPE, which is an input to the NPF package. ICELLTYPE == 0: *Confined cell* - Constant transmissivity ICELLTYPE != 0: *Convertible cell* - Transmissivity varies depending on the calculated head in the cell (based on the saturated cell thickness) In this example, all layers have a ICELLTYPE equal to 0, indicating confined cells. This is defined in the following section. Adding information to the NPF package -------------------------------------- The information for the NPF package is added to the gwf_model variable using `imod.mf6.NodePropertyFlow `_. The information included is the icelltype value (equal to zero), the array for the hydraulic conductivity (considered to be the same for the horizontal and vertical direction) and, optionally, the variable_vertical_conductance, dewatered, perched and save_flows options have been activated. For more details about the meaning of these variables and other variables available to be used within this package, please refer to the `documentation `_. .. GENERATED FROM PYTHON SOURCE LINES 119-130 .. code-block:: default gwf_model["npf"] = imod.mf6.NodePropertyFlow( icelltype=0, k=k, k33=k, variable_vertical_conductance=True, dewatered=True, perched=True, save_flows=True, ) .. GENERATED FROM PYTHON SOURCE LINES 131-145 Initial conditions package - IC ================================ This package reads the starting heads for a simulation. Starting heads interpolation ---------------------------- The starting heads to be used in this model are based on the interpolation of x-y head measurements, which were interpolated on a larger area. This example was created in this example --insert-link-here-- The heads were interpolated on a larger area, therefore these have to be clipped first .. GENERATED FROM PYTHON SOURCE LINES 145-162 .. code-block:: default initial = imod.data.hondsrug_initial() interpolated_head_larger = initial["head"] xmin = 237_500.0 xmax = 250_000.0 ymin = 559_000.0 ymax = 564_000.0 interpolated_head = interpolated_head_larger.sel( x=slice(xmin, xmax), y=slice(ymax, ymin) ) # Plotting the clipped interpolation fig, ax = plt.subplots() interpolated_head.plot.imshow(ax=ax) .. image-sg:: /examples/mf6/images/sphx_glr_hondsrug_002.png :alt: dx = 25.0, dy = -25.0 :srcset: /examples/mf6/images/sphx_glr_hondsrug_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 163-172 The final step is to assign the 2D heads interpolation to all the model layers (as a reference value) by using the xarray tool `xarray.full_like `_. The 3d idomain array is used as reference for the geometry and then its original values are replaced by NaNs. This array is combined with the interpolated_head array using the xarray `DataArray.combine_first `_ option. The final result is an starting_heads xarray where all layers have the 2d interpolated_head information. .. GENERATED FROM PYTHON SOURCE LINES 172-181 .. code-block:: default # Assign interpolated head values to all the model layers like_3d = xr.full_like(idomain, np.nan, dtype=float) starting_head = like_3d.combine_first(interpolated_head) # Consequently ensure no data is specified in inactive cells: starting_head = starting_head.where(idomain == 1) starting_head .. raw:: html 
<xarray.DataArray 'idomain' (layer: 13, y: 200, x: 500)>
    array([[[       nan,        nan,        nan, ..., 1.27008811,
             1.27823924, 1.28676046],
            [       nan,        nan,        nan, ..., 1.26336139,
             1.27183901, 1.28070347],
            [       nan,        nan,        nan, ..., 1.25620192,
             1.26504986, 1.27431479],
            ...,
            [7.28608515, 7.26550211, 7.24273429, ...,        nan,
                    nan,        nan],
            [7.2907364 , 7.26985601, 7.24672249, ...,        nan,
                    nan,        nan],
            [7.29728577, 7.27649525, 7.25355321, ...,        nan,
                    nan,        nan]],

           [[       nan,        nan,        nan, ...,        nan,
                    nan,        nan],
            [       nan,        nan,        nan, ...,        nan,
                    nan,        nan],
            [       nan,        nan,        nan, ...,        nan,
                    nan,        nan],
    ...
            [       nan,        nan,        nan, ...,        nan,
                    nan,        nan],
            [       nan,        nan,        nan, ...,        nan,
                    nan,        nan],
            [       nan,        nan,        nan, ...,        nan,
                    nan,        nan]],

           [[5.6706604 , 5.65319374, 5.63573351, ..., 1.27008811,
             1.27823924, 1.28676046],
            [5.67764711, 5.66021213, 5.64275829, ..., 1.26336139,
             1.27183901, 1.28070347],
            [5.68462064, 5.66720984, 5.64976361, ..., 1.25620192,
             1.26504986, 1.27431479],
            ...,
            [7.28608515, 7.26550211, 7.24273429, ..., 4.67051889,
             4.66997651, 4.66944365],
            [7.2907364 , 7.26985601, 7.24672249, ..., 4.70303061,
             4.70234778, 4.70168152],
            [7.29728577, 7.27649525, 7.25355321, ..., 4.73551559,
             4.73469247, 4.73389453]]])
    Coordinates:
      * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
      * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
        dx       float64 25.0
        dy       float64 -25.0
      * layer    (layer) int32 1 2 3 4 5 6 7 8 9 10 11 12 13


.. GENERATED FROM PYTHON SOURCE LINES 182-190 Adding information to the IC package ------------------------------------ The function for indicating the initial conditions is `imod.mf6.InitialConditions `_. It is necessary to indicate the value(s) to be considered as the initial (starting) head of the simulation. In this case, this value is equal to the previously created starting_head array. .. GENERATED FROM PYTHON SOURCE LINES 190-193 .. code-block:: default gwf_model["ic"] = imod.mf6.InitialConditions(starting_head) .. GENERATED FROM PYTHON SOURCE LINES 194-207 Constant head package - CHD =========================== This package allows to indicate if the head varies with time, if it is constant or if it is inactive. Constant head edge ------------------- The previously interpolated starting_head array will be used to define the constant head value which will be used along the model boundaries. A function is defined to indicate the location of the outer edge (returning a boolean array). .. GENERATED FROM PYTHON SOURCE LINES 207-216 .. code-block:: default def outer_edge(da): data = da.copy() from_edge = scipy.ndimage.binary_erosion(data) is_edge = (data == 1) & (from_edge == 0) return is_edge.astype(bool) .. GENERATED FROM PYTHON SOURCE LINES 217-221 For the next calculations, it is necessary to create a template array which can be used for assigning the corresponding geometry to other arrays. In this case, a 2d template is created based on the idomain layer information and filled with ones. .. GENERATED FROM PYTHON SOURCE LINES 221-225 .. code-block:: default like_2d = xr.full_like(idomain.isel(layer=0), 1) like_2d .. raw:: html 
<xarray.DataArray 'idomain' (y: 200, x: 500)>
    array([[1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           ...,
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1]])
    Coordinates:
      * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
      * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
        dx       float64 ...
        dy       float64 ...
        layer    int32 1


.. GENERATED FROM PYTHON SOURCE LINES 226-228 Using the previously created function and the 2d template, the outer edge is defined for this example. .. GENERATED FROM PYTHON SOURCE LINES 228-231 .. code-block:: default edge = outer_edge(xr.full_like(like_2d.drop_vars("layer"), 1)) .. GENERATED FROM PYTHON SOURCE LINES 232-249 Adding information to the CHD package -------------------------------------- To add the information to the CHD package within the gwf_model variable, the `imod.mf6.ConstantHead `_ function is used. The required information is the head array for this boundary condition. In this example, the starting_head array is selected where the idomain is > 0 (active) and it is located in the edge array. It is also possible (and optional) to indicate if the CHD information will be written to the listing file after it is read (print_input), if the constant head flow rates will be printed to the listing file for every stress period time step in which “BUDGET PRINT” is specified in Output Control (print_flows) and if the constant head flow terms will be written to the file specified with “BUDGET FILEOUT” in Output Control (save_flows). By default, these three options are set to False. .. GENERATED FROM PYTHON SOURCE LINES 249-257 .. code-block:: default gwf_model["chd"] = imod.mf6.ConstantHead( starting_head.where((idomain > 0) & edge), print_input=False, print_flows=True, save_flows=True, ) .. GENERATED FROM PYTHON SOURCE LINES 258-271 Recharge ======== This package is used to represent areally distributed recharge to the groundwater system. To calculate the recharge, the precipitation and evapotranspiration information from the KNMI website has been downloaded for the study area. This information is saved in netCDF files, which have been imported using the xarray function `xr.open_dataset `_, slicing the area to the model's miminum and maximum dimensions. Note that the meteorological data has mm/d as unit and this has to be converted to m/d for Modflow 6. .. GENERATED FROM PYTHON SOURCE LINES 272-285 .. code-block:: default xmin = 230_000.0 xmax = 257_000.0 ymin = 550_000.0 ymax = 567_000.0 meteorology = imod.data.hondsrug_meteorology() pp = meteorology["precipitation"] evt = meteorology["evapotranspiration"] pp = pp.sel(x=slice(xmin, xmax), y=slice(ymax, ymin)) / 1000.0 # from mm/d to m/d evt = evt.sel(x=slice(xmin, xmax), y=slice(ymax, ymin)) / 1000.0 # from mm/d to m/d .. GENERATED FROM PYTHON SOURCE LINES 286-296 Recharge - Steady state ----------------------- For the steady state conditions of the model, the data from the period 2000 to 2009 was considered as reference. The initial information was sliced to this time period and averaged to obtain the a mean value grid. This process was done for both precipitation and evapotranspiration datasets. **Precipitation** .. GENERATED FROM PYTHON SOURCE LINES 296-302 .. code-block:: default pp_ss = pp.sel(time=slice("2000-01-01", "2009-12-31")) pp_ss_mean = pp_ss.mean(dim="time") fig, ax = plt.subplots() pp_ss_mean.plot(ax=ax) .. image-sg:: /examples/mf6/images/sphx_glr_hondsrug_003.png :alt: dx = 1e+03, dy = -1e+03 :srcset: /examples/mf6/images/sphx_glr_hondsrug_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 303-304 **Evapotranspiration** .. GENERATED FROM PYTHON SOURCE LINES 304-310 .. code-block:: default evt_ss = evt.sel(time=slice("2000-01-01", "2009-12-31")) evt_ss_mean = evt_ss.mean(dim="time") fig, ax = plt.subplots() evt_ss_mean.plot(ax=ax) .. image-sg:: /examples/mf6/images/sphx_glr_hondsrug_004.png :alt: dx = 1e+03, dy = -1e+03 :srcset: /examples/mf6/images/sphx_glr_hondsrug_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 311-313 For the recharge calculation, a first estimate is the difference between the precipitation and evapotranspiration values. .. GENERATED FROM PYTHON SOURCE LINES 313-319 .. code-block:: default rch_ss = pp_ss_mean - evt_ss_mean fig, ax = plt.subplots() rch_ss.plot.imshow(ax=ax) .. image-sg:: /examples/mf6/images/sphx_glr_hondsrug_005.png :alt: dx = 1e+03, dy = -1e+03 :srcset: /examples/mf6/images/sphx_glr_hondsrug_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 320-325 Recharge - Transient -------------------- The transient model will encompass the period from 2010 to 2015. The initial pp and evt datasets have been sliced to this time frame. .. GENERATED FROM PYTHON SOURCE LINES 325-329 .. code-block:: default pp_trans = pp.sel(time=slice("2010-01-01", "2015-12-31")) evt_trans = evt.sel(time=slice("2010-01-01", "2015-12-31")) .. GENERATED FROM PYTHON SOURCE LINES 330-334 As previously done, it is assumed that the recharge is equal to the difference between precipitation and evapotranspiration as a first estimate. Furthermore, the negative values found after doing this calculation have been replaced by zeros, as the recharge should not have a negative value. .. GENERATED FROM PYTHON SOURCE LINES 334-338 .. code-block:: default rch_trans = pp_trans - evt_trans rch_trans = rch_trans.where(rch_trans > 0, 0) # check negative values .. GENERATED FROM PYTHON SOURCE LINES 339-342 The original information is on a daily step, so it is going to be resampled to a yearly step by using the xarray function `Dataset.resample `_. .. GENERATED FROM PYTHON SOURCE LINES 342-346 .. code-block:: default rch_trans_yr = rch_trans.resample(time="A", label="left").mean() rch_trans_yr .. raw:: html 
<xarray.DataArray (time: 6, y: 17, x: 27)>
    array([[[0.00191725, 0.00191589, 0.00191248, ..., 0.00158879,
             0.00160275, 0.00161772],
            [0.00193096, 0.00193016, 0.0019277 , ..., 0.00159777,
             0.00160856, 0.00162071],
            [0.00194592, 0.00194606, 0.00194384, ..., 0.00160863,
             0.00161658, 0.0016263 ],
            ...,
            [0.00205102, 0.0020491 , 0.00204451, ..., 0.00175834,
             0.00175477, 0.00175247],
            [0.0020481 , 0.00204493, 0.00203995, ..., 0.00177083,
             0.0017666 , 0.00176362],
            [0.00204444, 0.002041  , 0.00203658, ..., 0.00178251,
             0.00177807, 0.00177411]],

           [[0.00175232, 0.0017549 , 0.00175632, ..., 0.00160428,
             0.0016166 , 0.00163142],
            [0.00175832, 0.00176321, 0.00176716, ..., 0.00161909,
             0.00162744, 0.00163787],
            [0.00176593, 0.00177297, 0.00177865, ..., 0.00163613,
             0.00164029, 0.00164618],
    ...
            [0.00190413, 0.0019    , 0.00189647, ..., 0.00168478,
             0.00167463, 0.00166721],
            [0.00191761, 0.00191336, 0.00190866, ..., 0.00169624,
             0.00168655, 0.0016792 ],
            [0.0019308 , 0.00192603, 0.0019212 , ..., 0.00170712,
             0.00169771, 0.00169032]],

           [[0.00226806, 0.00225167, 0.00223808, ..., 0.00221006,
             0.00219559, 0.00218172],
            [0.0022924 , 0.00227892, 0.00226807, ..., 0.00220936,
             0.00219364, 0.00217812],
            [0.00231699, 0.00230603, 0.00229702, ..., 0.00221065,
             0.0021932 , 0.00217596],
            ...,
            [0.00249076, 0.00248442, 0.00247622, ..., 0.00224485,
             0.00222806, 0.00221361],
            [0.00249192, 0.00248449, 0.00247561, ..., 0.00224585,
             0.00223096, 0.00221743],
            [0.00249155, 0.00248376, 0.00247603, ..., 0.00224736,
             0.00223372, 0.00222133]]], dtype=float32)
    Coordinates:
      * x        (x) float64 2.305e+05 2.315e+05 2.325e+05 ... 2.555e+05 2.565e+05
      * y        (y) float64 5.665e+05 5.655e+05 5.645e+05 ... 5.515e+05 5.505e+05
        dx       float64 1e+03
        dy       float64 -1e+03
      * time     (time) datetime64[ns] 2009-12-31 2010-12-31 ... 2014-12-31


.. GENERATED FROM PYTHON SOURCE LINES 347-356 To create the final recharge for the transient simulation, the steady state information needs to be concatenated to the transient recharge data. The steady state simulation will be run for one second. This is achieved by using the numpy `Timedelta function `_, first creating a time delta of 1 second, which is assigned to the steady state recharge information. This dataset is then concatenated using the xarray function `xarray.concat `_ to the transient information and indicating that the dimension to join is "time". .. GENERATED FROM PYTHON SOURCE LINES 356-367 .. code-block:: default starttime = "2009-12-31" # Add first steady-state timedelta = np.timedelta64(1, "s") # 1 second duration for initial steady-state starttime_steady = np.datetime64(starttime) - timedelta rch_ss = rch_ss.assign_coords(time=starttime_steady) rch_ss_trans = xr.concat([rch_ss, rch_trans_yr], dim="time") rch_ss_trans .. raw:: html 
<xarray.DataArray (y: 17, x: 27, time: 7)>
    array([[[0.00091255, 0.00191725, 0.00175232, ..., 0.00167949,
             0.00161827, 0.00226806],
            [0.0009093 , 0.00191589, 0.0017549 , ..., 0.001689  ,
             0.00161518, 0.00225167],
            [0.00090591, 0.00191248, 0.00175632, ..., 0.00169961,
             0.00161301, 0.00223808],
            ...,
            [0.00078183, 0.00158879, 0.00160428, ..., 0.00174276,
             0.00162122, 0.00221006],
            [0.00078091, 0.00160275, 0.0016166 , ..., 0.00174685,
             0.0016207 , 0.00219559],
            [0.00078071, 0.00161772, 0.00163142, ..., 0.00175276,
             0.00162407, 0.00218172]],

           [[0.000913  , 0.00193096, 0.00175832, ..., 0.00168658,
             0.00164151, 0.0022924 ],
            [0.0009098 , 0.00193016, 0.00176321, ..., 0.00169678,
             0.00163919, 0.00227892],
            [0.0009065 , 0.0019277 , 0.00176716, ..., 0.00170842,
             0.00163841, 0.00226807],
    ...
            [0.00085378, 0.00177083, 0.00186831, ..., 0.00174713,
             0.00169624, 0.00224585],
            [0.00084055, 0.0017666 , 0.00185661, ..., 0.00173329,
             0.00168655, 0.00223096],
            [0.00082746, 0.00176362, 0.00184691, ..., 0.00172136,
             0.0016792 , 0.00221743]],

           [[0.00098196, 0.00204444, 0.00188612, ..., 0.00182747,
             0.0019308 , 0.00249155],
            [0.00098031, 0.002041  , 0.00189853, ..., 0.00183236,
             0.00192603, 0.00248376],
            [0.00097932, 0.00203658, 0.00191231, ..., 0.00183771,
             0.0019212 , 0.00247603],
            ...,
            [0.00086102, 0.00178251, 0.00188156, ..., 0.00175411,
             0.00170712, 0.00224736],
            [0.00084771, 0.00177807, 0.00187044, ..., 0.0017413 ,
             0.00169771, 0.00223372],
            [0.0008352 , 0.00177411, 0.00186046, ..., 0.00173005,
             0.00169032, 0.00222133]]], dtype=float32)
    Coordinates:
      * x        (x) float64 2.305e+05 2.315e+05 2.325e+05 ... 2.555e+05 2.565e+05
      * y        (y) float64 5.665e+05 5.655e+05 5.645e+05 ... 5.515e+05 5.505e+05
        dx       float64 1e+03
        dy       float64 -1e+03
      * time     (time) datetime64[ns] 2009-12-30T23:59:59 2009-12-31 ... 2014-12-31


.. GENERATED FROM PYTHON SOURCE LINES 368-377 The data obtained from KNMI has different grid dimensions than the one considered in this example. To fix this, imod-python includes the option `imod.prepare.Regridder `_, which modifies the original grid dimensions to a different one. It is also possible to define the regridding method such as ``nearest``, ``multilinear``, ``mean``, among others. In this case, ``mean`` was selected and the 2d template (like_2d) was used as reference, as this is the geometry to be considered in the model. .. GENERATED FROM PYTHON SOURCE LINES 377-381 .. code-block:: default rch_ss_trans = imod.prepare.Regridder(method="mean").regrid(rch_ss_trans, like_2d) rch_ss_trans .. raw:: html 
<xarray.DataArray (y: 200, x: 500, time: 7)>
    array([[[0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
             0.00169166, 0.00232221],
            [0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
             0.00169166, 0.00232221],
            [0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
             0.00169166, 0.00232221],
            ...,
            [0.00083373, 0.00163785, 0.00171206, ..., 0.00176909,
             0.00164895, 0.00231468],
            [0.00083373, 0.00163785, 0.00171206, ..., 0.00176909,
             0.00164895, 0.00231468],
            [0.00083373, 0.00163785, 0.00171206, ..., 0.00176909,
             0.00164895, 0.00231468]],

           [[0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
             0.00169166, 0.00232221],
            [0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
             0.00169166, 0.00232221],
            [0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
             0.00169166, 0.00232221],
    ...
            [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
             0.00171221, 0.0023817 ],
            [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
             0.00171221, 0.0023817 ],
            [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
             0.00171221, 0.0023817 ]],

           [[0.00091826, 0.00198224, 0.00191962, ..., 0.00183427,
             0.00179237, 0.00242177],
            [0.00091826, 0.00198224, 0.00191962, ..., 0.00183427,
             0.00179237, 0.00242177],
            [0.00091826, 0.00198224, 0.00191962, ..., 0.00183427,
             0.00179237, 0.00242177],
            ...,
            [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
             0.00171221, 0.0023817 ],
            [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
             0.00171221, 0.0023817 ],
            [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
             0.00171221, 0.0023817 ]]])
    Coordinates:
      * time     (time) datetime64[ns] 2009-12-30T23:59:59 2009-12-31 ... 2014-12-31
      * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
      * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
        dx       float64 25.0
        dy       float64 -25.0
        layer    int32 1


.. GENERATED FROM PYTHON SOURCE LINES 382-388 The previously created recharge array is a 2D array that needs to be assigned to a 3D array. This is done using the xarray `DataArray.where `_ option, where the recharge values are applied to the cells where the idomain value is larger than zero (that is, the active cells) and for the uppermost active cell (indicated by the minimum layer number). .. GENERATED FROM PYTHON SOURCE LINES 388-394 .. code-block:: default rch_total = rch_ss_trans.where( idomain["layer"] == idomain["layer"].where(idomain > 0).min("layer") ) rch_total .. raw:: html 
<xarray.DataArray (y: 200, x: 500, time: 7, layer: 13)>
    array([[[[       nan,        nan, 0.00089943, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.00190719, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.00182535, ...,        nan,
                     nan,        nan],
             ...,
             [       nan,        nan, 0.00179266, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.00169166, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.00232221, ...,        nan,
                     nan,        nan]],

            [[       nan,        nan, 0.00089943, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.00190719, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.00182535, ...,        nan,
                     nan,        nan],
    ...
             [       nan,        nan, 0.00180298, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.00171221, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.0023817 , ...,        nan,
                     nan,        nan]],

            [[       nan,        nan, 0.00088318, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.00169605, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.00183059, ...,        nan,
                     nan,        nan],
             ...,
             [       nan,        nan, 0.00180298, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.00171221, ...,        nan,
                     nan,        nan],
             [       nan,        nan, 0.0023817 , ...,        nan,
                     nan,        nan]]]])
    Coordinates:
      * time     (time) datetime64[ns] 2009-12-30T23:59:59 2009-12-31 ... 2014-12-31
      * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
      * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
        dx       float64 25.0
        dy       float64 -25.0
      * layer    (layer) int32 1 2 3 4 5 6 7 8 9 10 11 12 13


.. GENERATED FROM PYTHON SOURCE LINES 395-398 Finally, transposing the array dimensions using `DataArray.transpose `_ so they are in the correct order. .. GENERATED FROM PYTHON SOURCE LINES 398-405 .. code-block:: default rch_total = rch_total.transpose("time", "layer", "y", "x") rch_total fig, ax = plt.subplots() rch_total.isel(layer=2, time=6).plot.imshow(ax=ax) .. image-sg:: /examples/mf6/images/sphx_glr_hondsrug_006.png :alt: time = 2014-12-31, dx = 25.0, dy = -25.0, layer... :srcset: /examples/mf6/images/sphx_glr_hondsrug_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 406-413 Adding information to the RCH package -------------------------------------- The information for the RCH package is added with the function `imod.mf6.Recharge `_. It is required to insert the recharge flux rate, and it is optional to include the print_input, print_flows and save_flows information. .. GENERATED FROM PYTHON SOURCE LINES 413-416 .. code-block:: default gwf_model["rch"] = imod.mf6.Recharge(rch_total) .. GENERATED FROM PYTHON SOURCE LINES 417-429 Drainage package - DRN ======================= The drain package is used to simulate features that remove water from the aquifer, such as agricultural drains or springs. This occurs at a rate proportional to the head difference between the head in the aquifer and the drain elevation (the head in the aquifer has to be above that elevation). The conductance is the proportionality constant. Import drainage information ---------------------------- .. GENERATED FROM PYTHON SOURCE LINES 429-437 .. code-block:: default drainage = imod.data.hondsrug_drainage() pipe_cond = drainage["conductance"] pipe_elev = drainage["elevation"] pipe_cond .. raw:: html 
<xarray.DataArray 'conductance' (layer: 4, y: 200, x: 500)>
    [400000 values with dtype=float32]
    Coordinates:
      * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
      * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
        dx       float64 ...
        dy       float64 ...
      * layer    (layer) int32 1 3 5 6


.. GENERATED FROM PYTHON SOURCE LINES 438-448 Adding information to the DRN package ------------------------------------- To add the information to the DRN package within the gwf_model variable, the `imod.mf6.Drainage `_ function is used. It is required to add the previously created arrays for the drain elevation and the drain conductance. It is optional to insert the information for ``print_input``, ``print_flows`` and ``save_flows`` which are set to False by default. .. GENERATED FROM PYTHON SOURCE LINES 448-451 .. code-block:: default gwf_model["drn-pipe"] = imod.mf6.Drainage(conductance=pipe_cond, elevation=pipe_elev) .. GENERATED FROM PYTHON SOURCE LINES 452-460 River package - RIV =================== This package simulates the effects of flow between surface-water features and groundwater systems. Import river information ------------------------ .. GENERATED FROM PYTHON SOURCE LINES 460-466 .. code-block:: default river = imod.data.hondsrug_river() riv_cond = river["conductance"] riv_stage = river["stage"] riv_bot = river["bottom"] .. GENERATED FROM PYTHON SOURCE LINES 467-473 Adding information to the RIV package ------------------------------------- The data is assigned to the gwf_model variable by using `imod.mf6.River `_, based on the previously imported conductance, stage and bottom arrays. .. GENERATED FROM PYTHON SOURCE LINES 473-478 .. code-block:: default gwf_model["riv"] = imod.mf6.River( conductance=riv_cond, stage=riv_stage, bottom_elevation=riv_bot ) .. GENERATED FROM PYTHON SOURCE LINES 479-484 Storage package - STO ====================== When the STO Package is included in a model, storage changes will be calculated, and thus, the model will be transient. .. GENERATED FROM PYTHON SOURCE LINES 484-509 .. code-block:: default ss = 0.0003 layer = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]) sy = xr.DataArray( [0.16, 0.16, 0.16, 0.16, 0.15, 0.15, 0.15, 0.15, 0.14, 0.14, 0.14, 0.14, 0.14], {"layer": layer}, ("layer",), ) times_sto = np.array( [ "2009-12-30T23:59:59.00", "2009-12-31T00:00:00.00", "2010-12-31T00:00:00.00", "2011-12-31T00:00:00.00", "2012-12-31T00:00:00.00", "2013-12-31T00:00:00.00", "2014-12-31T00:00:00.00", ], dtype="datetime64[ns]", ) transient = xr.DataArray( [False, True, True, True, True, True, True], {"time": times_sto}, ("time",) ) .. GENERATED FROM PYTHON SOURCE LINES 510-517 Adding information to the STO package ------------------------------------- The data is assigned to the gwf_model variable by using `imod.mf6.Storage `_. It is necessary to indicate the values of specific storage, specific yield and if the layers are convertible. .. GENERATED FROM PYTHON SOURCE LINES 517-526 .. code-block:: default gwf_model["sto"] = imod.mf6.SpecificStorage( specific_storage=ss, specific_yield=sy, transient=transient, convertible=0, save_flows=True, ) .. GENERATED FROM PYTHON SOURCE LINES 527-543 Output Control package - OC =========================== This package determines how and when heads are printed to the listing file and/or written to a separate binary output file Adding information to the OC package ------------------------------------ The function `imod.mf6.OutputControl `_ is used to store the information for this package. It is possible to indicate if the heads and budget information is saved at the end of each stress period (``last``), for all timesteps a stress period (``all``), or at the start of a stress period (``first``) .. GENERATED FROM PYTHON SOURCE LINES 543-546 .. code-block:: default gwf_model["oc"] = imod.mf6.OutputControl(save_head="last", save_budget="last") .. GENERATED FROM PYTHON SOURCE LINES 547-559 Model simulation ================ In MODFLOW 6, the concept of "model" is that part of the program that solves a hydrologic process. MODFLOW 6 documentation supports one type of model: the GWF Model. It is possible within the MODFLOW 6 framewotk to solve multiple, tightly coupled, numerical models in a single system of equation, which may be multiple models of the same type or of different types. The previously created gwf_model variable now contains the information from all the variables. .. GENERATED FROM PYTHON SOURCE LINES 559-562 .. code-block:: default gwf_model .. rst-class:: sphx-glr-script-out .. code-block:: none {'dis': , 'npf': , 'ic': , 'chd': , 'rch': , 'drn-pipe': , 'riv': , 'sto': , 'oc': } .. GENERATED FROM PYTHON SOURCE LINES 563-569 Attach the model information to a simulation -------------------------------------------- The function `imod.mf6.Modflow6Simulation `_ allows to assign models to a simulation (in this case, the gwf_model). .. GENERATED FROM PYTHON SOURCE LINES 569-573 .. code-block:: default simulation = imod.mf6.Modflow6Simulation("mf6-mipwa2-example") simulation["GWF_1"] = gwf_model .. GENERATED FROM PYTHON SOURCE LINES 574-580 Solver settings --------------- The solver settings are indicated using `imod.mf6.Solution `_. If the values are not indicated manually, the defaults values will be considered. .. GENERATED FROM PYTHON SOURCE LINES 580-598 .. code-block:: default simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", csv_output=False, no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, inner_dvclose=1.0e-4, inner_rclose=0.001, inner_maximum=100, linear_acceleration="cg", scaling_method=None, reordering_method=None, relaxation_factor=0.97, ) .. GENERATED FROM PYTHON SOURCE LINES 599-603 Assign time discretization -------------------------- The time discretization of this model is 6 years. .. GENERATED FROM PYTHON SOURCE LINES 603-608 .. code-block:: default simulation.create_time_discretization( additional_times=["2009-12-30T23:59:59.000000000", "2015-12-31T00:00:00.000000000"] ) .. GENERATED FROM PYTHON SOURCE LINES 609-617 Run the model ------------- .. note:: The following lines assume the ``mf6`` executable is available on your PATH. :ref:`The Modflow 6 examples introduction ` shortly describes how to add it to yours. .. GENERATED FROM PYTHON SOURCE LINES 617-622 .. code-block:: default modeldir = imod.util.temporary_directory() simulation.write(modeldir, binary=False, validate=False) simulation.run() .. GENERATED FROM PYTHON SOURCE LINES 623-634 Results visualization ===================== The next section indicated how to visualize the model results. Import heads results -------------------- The heads results are imported using `imod.mf6.open_hds `_ and indicating the location of the heads, as well as the DIS file. .. GENERATED FROM PYTHON SOURCE LINES 634-640 .. code-block:: default hds = imod.mf6.open_hds( modeldir / "GWF_1" / "GWF_1.hds", modeldir / "GWF_1" / "dis.dis.grb" ) hds .. raw:: html 
<xarray.DataArray 'head' (time: 7, layer: 13, y: 200, x: 500)>
    dask.array<stack, shape=(7, 13, 200, 500), dtype=float64, chunksize=(1, 13, 200, 500), chunktype=numpy.ndarray>
    Coordinates:
      * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
      * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
        dx       float64 25.0
        dy       float64 -25.0
      * layer    (layer) int64 1 2 3 4 5 6 7 8 9 10 11 12 13
      * time     (time) float64 1.157e-05 365.0 730.0 ... 1.826e+03 2.191e+03


.. GENERATED FROM PYTHON SOURCE LINES 641-642 We can plot the data of an individual layer as follows .. GENERATED FROM PYTHON SOURCE LINES 642-645 .. code-block:: default fig, ax = plt.subplots() hds.sel(layer=3).isel(time=3).plot(ax=ax) .. image-sg:: /examples/mf6/images/sphx_glr_hondsrug_007.png :alt: dx = 25.0, dy = -25.0, layer = 3, time = 1.096e+03 :srcset: /examples/mf6/images/sphx_glr_hondsrug_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 646-649 As you can see layer 3 has some missing cells in the west Whereas layer 4 only contains active cells in the eastern peatland area .. GENERATED FROM PYTHON SOURCE LINES 649-651 .. code-block:: default fig, ax = plt.subplots() hds.sel(layer=4).isel(time=3).plot(ax=ax) .. image-sg:: /examples/mf6/images/sphx_glr_hondsrug_008.png :alt: dx = 25.0, dy = -25.0, layer = 4, time = 1.096e+03 :srcset: /examples/mf6/images/sphx_glr_hondsrug_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 652-654 Layer 5 contains more data towards the west, but has no active cells in the centre. .. GENERATED FROM PYTHON SOURCE LINES 654-657 .. code-block:: default fig, ax = plt.subplots() hds.sel(layer=5).isel(time=3).plot(ax=ax) .. image-sg:: /examples/mf6/images/sphx_glr_hondsrug_009.png :alt: dx = 25.0, dy = -25.0, layer = 5, time = 1.096e+03 :srcset: /examples/mf6/images/sphx_glr_hondsrug_009.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 658-669 As you can see the data is individual layers have lots of inactive in different places. It is difficult for this model to get a good idea what is happening across the area based on 1 layer alone. Luckily xarray allows us to compute the mean across a selection of layers and plot this. By first selecting 3 layers with ``sel```, and then computing the mean across the layer dimension with ``mean(dim="layer")``. .. GENERATED FROM PYTHON SOURCE LINES 669-673 .. code-block:: default fig, ax = plt.subplots() hds.sel(layer=slice(3, 5)).mean(dim="layer").isel(time=3).plot(ax=ax) .. image-sg:: /examples/mf6/images/sphx_glr_hondsrug_010.png :alt: dx = 25.0, dy = -25.0, time = 1.096e+03 :srcset: /examples/mf6/images/sphx_glr_hondsrug_010.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 674-680 Assign dates to head -------------------- MODFLOW6 has no concept of a calendar, so the output is not labelled only in terms of "time since start" in floating point numbers. For this model the time unit is days and we can assign a date coordinate as follows: .. GENERATED FROM PYTHON SOURCE LINES 680-685 .. code-block:: default starttime = pd.to_datetime("2000-01-01") timedelta = pd.to_timedelta(hds["time"], "D") hds = hds.assign_coords(time=starttime + timedelta) .. GENERATED FROM PYTHON SOURCE LINES 686-692 Extract head at points ---------------------- A typical operation is to extract simulated heads at point locations to compare them with measurements. In this example, we select the heads at two points: .. GENERATED FROM PYTHON SOURCE LINES 692-697 .. code-block:: default x = [240_000.0, 244_000.0] y = [560_000.0, 562_000.0] selection = imod.select.points_values(hds, x=x, y=y) .. GENERATED FROM PYTHON SOURCE LINES 698-700 The result can be converted into a pandas dataframe for timeseries analysis, or written to a variety of tabular file formats. .. GENERATED FROM PYTHON SOURCE LINES 700-704 .. code-block:: default dataframe = selection.to_dataframe().reset_index() dataframe = dataframe.rename(columns={"index": "id"}) dataframe .. raw:: html
time layer id x y dx dy head
0 2000-01-01 00:00:00.999993600 1 0 240012.5 560012.5 25.0 -25.0 NaN
1 2000-01-01 00:00:00.999993600 1 1 244012.5 562012.5 25.0 -25.0 10.534013
2 2000-01-01 00:00:00.999993600 2 0 240012.5 560012.5 25.0 -25.0 NaN
3 2000-01-01 00:00:00.999993600 2 1 244012.5 562012.5 25.0 -25.0 NaN
4 2000-01-01 00:00:00.999993600 3 0 240012.5 560012.5 25.0 -25.0 7.682111
... ... ... ... ... ... ... ... ...
177 2005-12-31 00:00:00.999993600 11 1 244012.5 562012.5 25.0 -25.0 4.787419
178 2005-12-31 00:00:00.999993600 12 0 240012.5 560012.5 25.0 -25.0 NaN
179 2005-12-31 00:00:00.999993600 12 1 244012.5 562012.5 25.0 -25.0 NaN
180 2005-12-31 00:00:00.999993600 13 0 240012.5 560012.5 25.0 -25.0 7.498579
181 2005-12-31 00:00:00.999993600 13 1 244012.5 562012.5 25.0 -25.0 4.787407

182 rows × 8 columns



.. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 11.814 seconds) .. _sphx_glr_download_examples_mf6_hondsrug.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: hondsrug.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: hondsrug.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_