From 1ab55d90ba49277db2f5071f44c00c31a55fcc4b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Wed, 27 Sep 2023 17:28:08 +0200 Subject: [PATCH 01/11] Add first version of uzf! --- nlmod/dims/grid.py | 25 ++++ nlmod/dims/time.py | 2 + nlmod/gwf/gwf.py | 25 ++++ nlmod/gwf/recharge.py | 341 ++++++++++++++++++++++++++++++++++++------ nlmod/read/knmi.py | 2 +- tests/test_020_uzf.py | 50 +++++++ 6 files changed, 396 insertions(+), 49 deletions(-) create mode 100644 tests/test_020_uzf.py diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index d95d0e64..049e8a4b 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -840,6 +840,31 @@ def lcid_to_reclist( return reclist +def cols_to_reclist(ds, cellids, *args, cellid_column=0): + """Create a reclist for stress period data from a set of cellids. + + Parameters + ---------- + ds : xarray.Dataset + dataset with model data. Should have dimensions (layer, icell2d). + cellids : tuple of length 2 or 3 + tuple with indices of the cells that will be used to create the list. For a + structured grid, cellids represents (layer, row, column). For a vertex grid + cellid reprsents (layer, icell2d). + args : xarray.DatArray, str, int or float + the args parameter represents the data to be used as the columns in the reclist. + See col_to_list of the allowed values. + cellid_column : int, optional + Adds the cellid ((layer, row, col) or (layer, icell2d)) to the reclist in this + column number. Do not add cellid when cellid_column is None. The default is 0. + + """ + cols = [col_to_list(col, ds, cellids) for col in args] + if cellid_column is not None: + cols.insert(cellid_column, list(zip(*cellids))) + return list(zip(*cols)) + + def da_to_reclist( ds, mask, diff --git a/nlmod/dims/time.py b/nlmod/dims/time.py index 22bc9024..697b5bbe 100644 --- a/nlmod/dims/time.py +++ b/nlmod/dims/time.py @@ -509,6 +509,8 @@ def dataframe_to_flopy_timeseries( timeseries = [(i,) + tuple(v) for i, v in zip(df.index, df.values)] if package is None: return timeseries + if filename is None: + filename = f"{package.filename}_ts" if time_series_namerecord is None: time_series_namerecord = list(df.columns) diff --git a/nlmod/gwf/gwf.py b/nlmod/gwf/gwf.py index 62d478ba..8653fbe0 100644 --- a/nlmod/gwf/gwf.py +++ b/nlmod/gwf/gwf.py @@ -854,6 +854,31 @@ def evt(ds, gwf, pname="evt", **kwargs): return evt +def uzf(ds, gwf, pname="uzf", **kwargs): + """create unsaturated zone flow package from model dataset. + + Parameters + ---------- + ds : xarray.Dataset + dataset with model data. + gwf : flopy ModflowGwf + groundwaterflow object. + pname : str, optional + package name + + Returns + ------- + uzf : flopy ModflowGwfuzf + uzf package + """ + logger.info("creating mf6 UZF") + + # create uzf package + evt = recharge.ds_to_uzf(gwf, ds, pname=pname, **kwargs) + + return evt + + def _set_record(out, budget, output="head"): record = [] if isinstance(out, bool): diff --git a/nlmod/gwf/recharge.py b/nlmod/gwf/recharge.py index 751751dc..e18e5435 100644 --- a/nlmod/gwf/recharge.py +++ b/nlmod/gwf/recharge.py @@ -2,15 +2,18 @@ import flopy import numpy as np -from tqdm import tqdm +import pandas as pd +import xarray as xr -from ..dims.grid import da_to_reclist -from ..sim.sim import get_tdis_perioddata +from ..dims.grid import da_to_reclist, cols_to_reclist +from ..dims.layers import get_idomain, get_first_active_layer_from_idomain +from ..dims.time import dataframe_to_flopy_timeseries +from ..util import _get_value_from_ds_datavar logger = logging.getLogger(__name__) -def ds_to_rch(gwf, ds, mask=None, pname="rch", **kwargs): +def ds_to_rch(gwf, ds, mask=None, pname="rch", recharge="recharge", **kwargs): """Convert the recharge data in the model dataset to a rch package with time series. @@ -27,7 +30,7 @@ def ds_to_rch(gwf, ds, mask=None, pname="rch", **kwargs): Returns ------- - rch : flopy.mf6.modflow.mfgwfrch.ModflowGwfrch + rch : flopy.mf6.ModflowGwfrch recharge package """ # check for nan values @@ -36,19 +39,18 @@ def ds_to_rch(gwf, ds, mask=None, pname="rch", **kwargs): # get stress period data if ds["steady"].all(): - recharge = "recharge" - if "time" in ds["recharge"].dims: - mask = ds["recharge"].isel(time=0) != 0 + if "time" in ds[recharge].dims: + mask = ds[recharge].isel(time=0) != 0 else: - mask = ds["recharge"] != 0 + mask = ds[recharge] != 0 else: + rch_name_arr, rch_unique_dic = _get_unique_series(ds, recharge, pname) recharge = "rch_name" - rch_name_arr, rch_unique_dic = _get_unique_series(ds, "recharge", pname) ds[recharge] = ds["top"].dims, rch_name_arr if mask is not None: - mask = (ds["rch_name"] != "") & mask + mask = (ds[recharge] != "") & mask else: - mask = ds["rch_name"] != "" + mask = ds[recharge] != "" spd = da_to_reclist( ds, @@ -78,7 +80,9 @@ def ds_to_rch(gwf, ds, mask=None, pname="rch", **kwargs): return rch -def ds_to_evt(gwf, ds, pname="evt", nseg=1, surface=None, depth=None, **kwargs): +def ds_to_evt( + gwf, ds, pname="evt", rate="evaporation", nseg=1, surface=None, depth=None, **kwargs +): """Convert the evaporation data in the model dataset to a evt package with time series. @@ -110,7 +114,7 @@ def ds_to_evt(gwf, ds, pname="evt", nseg=1, surface=None, depth=None, **kwargs): Returns ------- - evt : flopy.mf6.modflow.mfgwfevt.ModflowGwfevt + evt : flopy.mf6.ModflowGwfevt evapotranspiiration package. """ assert nseg == 1, "More than one evaporation segment not yet supported" @@ -124,18 +128,17 @@ def ds_to_evt(gwf, ds, pname="evt", nseg=1, surface=None, depth=None, **kwargs): depth = 1.0 # check for nan values - if ds["evaporation"].isnull().any(): + if ds[rate].isnull().any(): raise ValueError("please remove nan values in evaporation data array") # get stress period data if ds["steady"].all(): - if "time" in ds["evaporation"].dims: - mask = ds["evaporation"].isel(time=0) != 0 + if "time" in ds[rate].dims: + mask = ds[rate].isel(time=0) != 0 else: - mask = ds["evaporation"] != 0 - rate = "evaporation" + mask = ds[rate] != 0 else: - evt_name_arr, evt_unique_dic = _get_unique_series(ds, "evaporation", pname) + evt_name_arr, evt_unique_dic = _get_unique_series(ds, rate, pname) ds["evt_name"] = ds["top"].dims, evt_name_arr mask = ds["evt_name"] != "" @@ -172,6 +175,263 @@ def ds_to_evt(gwf, ds, pname="evt", nseg=1, surface=None, depth=None, **kwargs): return evt +def ds_to_uzf( + gwf, + ds, + mask=None, + pname="uzf", + surfdep=0.05, + vks="kv", + thtr=0.1, + thts=0.3, + thti=0.1, + eps=3.5, + landflag=None, + finf="recharge", + pet="evaporation", + extdp=None, + extwc=None, + ha=None, + hroot=None, + rootact=None, + simulate_et=True, + linear_gwet=True, + unsat_etwc=False, + unsat_etae=False, + **kwargs, +): + """Create a unsaturated zone flow package for modflow 6. This method adds uzf-cells + to all active Modflow cells (unless mask is specified). + + Parameters + ---------- + gwf : flopy.mf6.modflow.mfgwf.ModflowGwf + groundwater flow model. + ds : xr.DataSet + dataset containing relevant model grid information + mask : xr.DataArray + data array containing mask, uzf is only added where mask is True + pname : str, optional + package name. The default is 'uzf'. + surfdep : float, str or array-like + surfdep is the surface depression depth of the UZF cell. When passed as string, + the array is obtained from ds. The default is 0.05 m. + vks : float, str or array-like + vks is the saturated vertical hydraulic conductivity of the UZF cell. This value + is used with the Brooks-Corey function and the simulated water content to + calculate the partially saturated hydraulic conductivity. When passed as string, + the array is obtained from ds. The default is 'kv'. + thtr : float, str or array-like + thtr is the residual (irreducible) water content of the UZF cell. This residual + water is not available to plants and will not drain into underlying aquifer + cells. When passed as string, the array is obtained from ds. The default is + 0.1. + thts : float, str or array-like + thts is the saturated water content of the UZF cell. The values for saturated + and residual water content should be set in a manner that is consistent with the + specific yield value specified in the Storage Package. The saturated water + content must be greater than the residual content. When passed as string, the + array is obtained from ds. The default is 0.3. + thti : float, str or array-like + thti is the initial water content of the UZF cell. The value must be greater + than or equal to the residual water content and less than or equal to the + saturated water content. When passed as string, the array is obtained from ds. + The default is 0.1. + eps : float, str or array-like + eps is the exponent used in the Brooks-Corey function. The Brooks-Corey function + is used by UZF to calculated hydraulic conductivity under partially saturated + conditions as a function of water content and the user-specified saturated + hydraulic conductivity. Values must be between 3.5 and 14.0. When passed as + string, the array is obtained from ds. The default is 3.5. + landflag : xr.DataArray, optional + A DataArray with integer values, set to one for land surface cells indicating + that boundary conditions can be applied and data can be specified in the PERIOD + block. A value of 0 specifies a non-land surface cell. Landflag is determined + from ds when it is None. The default is None. + finf :float, str or array-like + The applied infiltration rate of the UZF cell (:math:`LT^{-1}`). When passed as + string, the array is obtained from ds. The default is "recharge". + pet : float, str or array-like + The potential evapotranspiration rate of the UZF cell and specified GWF cell. + Evapotranspiration is first removed from the unsaturated zone and any remaining + potential evapotranspiration is applied to the saturated zone. only + used if SIMULATE_ET. When passed as string, the array is obtained from ds. The + default is "evaporation". + extdp : float, optional + Value that defines the evapotranspiration extinction depth of the UZF cell, in + meters below the top of the model. Set to 2.0 meter when None. The default is + None. + extwc : float, optional + The evapotranspiration extinction water content of the UZF cell. Only used if + SIMULATE_ET and UNSAT_ETWC. Set to thtr when None. The default is None. + ha : float, optional + The air entry potential (head) of the UZF cell. Only used if SIMULATE_ET and + UNSAT_ETAE. Set to 0.0 when None. The default is None. + hroot : float, optional + The root potential (head) of the UZF cell. Only used if SIMULATE_ET and + UNSAT_ETAE. Set to 0.0 when None. The default is None. + rootact : float, optional + the root activity function of the UZF cell. ROOTACT is the length of roots in + a given volume of soil divided by that volume. Values range from 0 to about 3 + :math:`cm^{-2}`, depending on the plant community and its stage of development. + Only used if SIMULATE_ET and UNSAT_ETAE. Set to 0.0 when None. The default is + None. + simulate_et : bool, optional + If True, ET in the unsaturated (UZF) and saturated zones (GWF) will be + simulated. The default is True. + linear_gwet : bool, optional + If True, groundwater ET will be simulated using the original ET formulation of + MODFLOW-2005. When False, and square_gwet is True as an extra argument, no + evaporation from the saturated zone will be simulated. The default is True. + unsat_etwc : bool, optional + If True, ET in the unsaturated zone will be simulated as a function of the + specified PET rate while the water content (THETA) is greater than the ET + extinction water content (EXTWC). The default is False. + unsat_etae : bool, optional + If True, ET in the unsaturated zone will be simulated using a capillary pressure + based formulation. Capillary pressure is calculated using the Brooks-Corey + retention function. The default is False. + ** kwargs : dict + Kwargs are passed onto flopy.mf6.ModflowGwfuzf + + + Returns + ------- + uzf : flopy.mf6.ModflowGwfuzf + Unsaturated zone flow package for Modflow 6. + """ + if mask is None: + mask = ds["area"] > 0 + + if "layer" not in mask.dims: + mask = mask.expand_dims(dim={"layer": ds.layer}) + + # only add uzf-cells in active cells + idomain = get_idomain(ds) + mask = mask & (idomain > 0) + + # generate packagedata + surfdep = _get_value_from_ds_datavar(ds, "surfdep", surfdep, return_da=True) + vks = _get_value_from_ds_datavar(ds, "vk", vks, return_da=True) + thtr = _get_value_from_ds_datavar(ds, "thtr", thtr, return_da=True) + thts = _get_value_from_ds_datavar(ds, "thts", thts, return_da=True) + thti = _get_value_from_ds_datavar(ds, "thti", thti, return_da=True) + eps = _get_value_from_ds_datavar(ds, "eps", eps, return_da=True) + + nuzfcells = int(mask.sum()) + cellids = np.where(mask) + + iuzno = xr.full_like(ds["botm"], -1, dtype=int) + iuzno.data[mask] = np.arange(nuzfcells) + + if landflag is None: + landflag = xr.full_like(ds["botm"], 0, dtype=int) + # set the landflag in the top layer to 1 + landflag[get_first_active_layer_from_idomain(idomain)] = 1 + + # determine ivertcon, by setting its value to iuzno of the layer below + ivertcon = xr.full_like(ds["botm"], -1, dtype=int) + ivertcon.data[:-1] = iuzno.data[1:] + # then use bfill to accont for inactive cells in the layer below, and set nans to -1 + ivertcon = ivertcon.where(ivertcon >= 0).bfill("layer").fillna(-1).astype(int) + + # packagedata : [iuzno, cellid, landflag, ivertcon, surfdep, vks, thtr, thts, thti, eps, boundname] + packagedata = cols_to_reclist( + ds, + cellids, + iuzno, + landflag, + ivertcon, + surfdep, + vks, + thtr, + thts, + thti, + eps, + cellid_column=1, + ) + + # add perioddata for all uzf cells that are at the surface + mask = landflag == 1 + + # perioddata : [iuzno, finf, pet, extdp, extwc, ha, hroot, rootact, aux] + finf_name_arr, uzf_unique_dic = _get_unique_series(ds, finf, "finf") + finf = "rch_name" + ds[finf] = ds["top"].dims, finf_name_arr + ds[finf] = ds[finf].expand_dims(dim={"layer": ds.layer}) + if mask is not None: + mask = (ds[finf] != "") & mask + else: + mask = ds[finf] != "" + + pet_name_arr, pet_unique_dic = _get_unique_series(ds, pet, "pet") + pet = "evt_name" + ds[pet] = ds["top"].dims, pet_name_arr + ds[pet] = ds[pet].expand_dims(dim={"layer": ds.layer}) + if mask is not None: + mask = (ds[pet] != "") & mask + else: + mask = ds[pet] != "" + + # combine the time series of finf and pet + uzf_unique_dic.update(pet_unique_dic) + + if extdp is None: + extdp = 2.0 + # EXTDP is always specified, but is only used if SIMULATE_ET + if simulate_et: + logger.info(f"Setting extinction depth (extdp) to {extdp} meter below top") + if extwc is None: + extwc = thtr + if simulate_et and unsat_etwc: + logger.info( + f"Setting evapotranspiration extinction water content (extwc) to {extwc}" + ) + if ha is None: + ha = 0.0 + if simulate_et and unsat_etae: + logger.info(f"Setting air entry potential (ha) to {ha}") + if hroot is None: + hroot = 0.0 + if simulate_et and unsat_etae: + logger.info(f"Setting root potential (hroot) to {hroot}") + if rootact is None: + rootact = 0.0 + if simulate_et and unsat_etae: + logger.info(f"Setting root activity function (rootact) to {rootact}") + + cellids = np.where(mask) + + perioddata = cols_to_reclist( + ds, + cellids, + iuzno, + finf, + pet, + extdp, + extwc, + ha, + hroot, + rootact, + cellid_column=None, + ) + + uzf = flopy.mf6.ModflowGwfuzf( + gwf, + nuzfcells=nuzfcells, + packagedata=packagedata, + perioddata={0: perioddata}, + simulate_et=simulate_et, + linear_gwet=linear_gwet, + unsat_etwc=unsat_etwc, + unsat_etae=unsat_etae, + **kwargs, + ) + + # create timeseries packages + _add_time_series(uzf, uzf_unique_dic, ds) + + def _get_unique_series(ds, var, pname): """Get the location and values of unique time series from a variable var in ds. @@ -247,31 +507,16 @@ def _add_time_series(package, rch_unique_dic, ds): ------- None. """ - # get timesteps - tdis_perioddata = get_tdis_perioddata(ds) - perlen_arr = [t[0] for t in tdis_perioddata] - time_steps_rch = [0.0] + np.array(perlen_arr).cumsum().tolist() - - for i, key in tqdm( - enumerate(rch_unique_dic.keys()), - total=len(rch_unique_dic.keys()), - desc="Building ts packages rch", - ): - # add extra time step to the time series object (otherwise flopy fails) - recharge_val = list(rch_unique_dic[key]) + [0.0] - - recharge = list(zip(time_steps_rch, recharge_val)) - if i == 0: - package.ts.initialize( - filename=f"{key}.ts", - timeseries=recharge, - time_series_namerecord=key, - interpolation_methodrecord="stepwise", - ) - else: - package.ts.append_package( - filename=f"{key}.ts", - timeseries=recharge, - time_series_namerecord=key, - interpolation_methodrecord="stepwise", - ) + # generate a DataFrame + df = pd.DataFrame(rch_unique_dic, index=ds.time) + if df.isna().any(axis=None): + # make sure there are no NaN's, as otherwise they will be filled by zeros later + raise (Exception("There cannot be nan's in the dictionary")) + # set the first value for the start-time as well + df.loc[pd.to_datetime(ds.time.start)] = df.iloc[0] + # combine the values with the start of each period, and ste the last value to 0.0 + df = df.sort_index().shift(-1).fillna(0.0) + + dataframe_to_flopy_timeseries( + df, ds=ds, package=package, interpolation_methodrecord="stepwise" + ) diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py index f104e819..92f530e0 100644 --- a/nlmod/read/knmi.py +++ b/nlmod/read/knmi.py @@ -38,7 +38,7 @@ def get_recharge(ds, method="linear", most_common_station=False): ---------- ds : xr.DataSet dataset containing relevant model grid information - method : bool, optional + method : str, optional If 'linear', calculate recharge by subtracting evaporation from precipitation. If 'separate', add precipitation as 'recharge' and evaporation as 'evaporation'. The defaults is 'linear'. diff --git a/tests/test_020_uzf.py b/tests/test_020_uzf.py new file mode 100644 index 00000000..03834201 --- /dev/null +++ b/tests/test_020_uzf.py @@ -0,0 +1,50 @@ +import numpy as np +import pandas as pd +import nlmod +import util + + +def test_uzf_structured(): + # %% create model Dataset + extent = [200_000, 202_000, 400_000, 403_000] + ds = util.get_ds_structured(extent, top=0, botm=np.linspace(-1, -10, 10)) + + time = pd.date_range("2022", "2023", freq="D") + ds = nlmod.time.set_ds_time(ds, start="2021", time=time) + + ds.update(nlmod.read.knmi.get_recharge(ds, method="separate")) + + # %% generate sim and gwf + # create simulation + sim = nlmod.sim.sim(ds) + + # create time discretisation + _ = nlmod.sim.tdis(ds, sim) + + # create groundwater flow model + gwf = nlmod.gwf.gwf(ds, sim) + + # create ims + _ = nlmod.sim.ims(sim) + + # Create discretization + _ = nlmod.gwf.dis(ds, gwf) + + # create node property flow + _ = nlmod.gwf.npf(ds, gwf) + + # Create the initial conditions package + _ = nlmod.gwf.ic(ds, gwf, starting_head=1.0) + + # Create the output control package + _ = nlmod.gwf.oc(ds, gwf) + + bhead = ds["botm"][1] + cond = ds["area"] * 1 + _ = nlmod.gwf.ghb(ds, gwf, bhead=bhead, cond=cond, layer=len(ds.layer) - 1) + + # create recharge package + _ = nlmod.gwf.uzf(ds, gwf) + + # %% run + _ = nlmod.sim.write_and_run(sim, ds) From f2c0c6ab1057f9de9bc1f347e06eabe3ee228fff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Fri, 29 Sep 2023 15:20:28 +0200 Subject: [PATCH 02/11] Add uzf observations and notebook --- docs/examples/17_unsaturated_zone_flow.ipynb | 415 +++++++++++++++++++ nlmod/gwf/recharge.py | 50 ++- tests/test_007_run_notebooks.py | 5 + tests/test_020_uzf.py | 2 +- 4 files changed, 466 insertions(+), 6 deletions(-) create mode 100644 docs/examples/17_unsaturated_zone_flow.ipynb diff --git a/docs/examples/17_unsaturated_zone_flow.ipynb b/docs/examples/17_unsaturated_zone_flow.ipynb new file mode 100644 index 00000000..57aa3c80 --- /dev/null +++ b/docs/examples/17_unsaturated_zone_flow.ipynb @@ -0,0 +1,415 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "71e9082c", + "metadata": {}, + "source": [ + "# Unsaturated zone flow\n", + "This notebook demonstrates the use of the Unsaturated Zone Flow (UZF) package in nlmod. The UZF-package can be used to simulate important processes in the unsaturated zone. These processes create a delay before precipitation reaches the saturated groundwater. In dry periods the water may even have evaporated before that. This notebook shows the difference in calculated head between a model with regular recharge (RCH) and evapotranspiration (EVT) packages, compared to a model with the UZF-package.\n", + "\n", + "We create a 1d model, of 1 row and 1 column, but with multiple layers, of a real location somewhere in the Netherlands. We use weather data from the KNMI as input for a transient simulation of 3 years with daily timetseps. This notebook can be used to vary the uzf-parameter, change the location (do not forget to alter the drain-elevation as well), or to play with the model-timestep." + ] + }, + { + "cell_type": "markdown", + "id": "0610edfd", + "metadata": {}, + "source": [ + "## Import packages and setup logging" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fbf2d8ce", + "metadata": {}, + "outputs": [], + "source": [ + "# import packages\n", + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import nlmod\n", + "\n", + "# set up pretty logging and show package versions\n", + "nlmod.util.get_color_logger()\n", + "nlmod.show_versions()" + ] + }, + { + "cell_type": "markdown", + "id": "9813a7da", + "metadata": {}, + "source": [ + "## Generate a model Dataset\n", + "We first set the model_name and model workspace, which we will use later to write the model files, and so we can determine the cache-directory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4cc92261", + "metadata": {}, + "outputs": [], + "source": [ + "model_name = \"test_uzf\"\n", + "model_ws = os.path.join(\"data\", model_name)\n", + "figdir, cachedir = nlmod.util.get_model_dirs(model_ws)" + ] + }, + { + "cell_type": "markdown", + "id": "8de98855", + "metadata": {}, + "source": [ + "We define a location with a corresponding drainage elevation. From the location we calculate an extent of 100 by 100 meter, and download REGIS." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "341f7f9d", + "metadata": {}, + "outputs": [], + "source": [ + "x = 37_400\n", + "y = 412_600\n", + "drn_elev = 4.0\n", + "\n", + "# round x and y to 100, so we will download only one cell in regis\n", + "x = np.floor(x / 100) * 100\n", + "y = np.floor(y / 100) * 100\n", + "dx = dy = 100\n", + "extent = [x, x + dx, y, y + dy]\n", + "regis = nlmod.read.regis.get_regis(extent, cachename='regis.nc', cachedir=cachedir)" + ] + }, + { + "cell_type": "markdown", + "id": "76deaeb7", + "metadata": {}, + "source": [ + "As the REGIS-data only contains one cell, we can visualize the properties of the layers in a pandas DataFrame. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bbd30a5a", + "metadata": {}, + "outputs": [], + "source": [ + "regis.sel(x=regis.x[0], y=regis.y[0]).to_pandas()" + ] + }, + { + "cell_type": "markdown", + "id": "cdb20bc9", + "metadata": {}, + "source": [ + "As you can see, there are some NaN-values in the hydaulic conductivities (`kh` and `kv`). These will be filled when making a model Dataset, usinge fill-values and anisotropy values. See the info-messages after the commands below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71466503", + "metadata": {}, + "outputs": [], + "source": [ + "ds = nlmod.to_model_ds(\n", + " regis,\n", + " model_name=model_name,\n", + " model_ws=model_ws,\n", + " fill_value_kh=10.0,\n", + " fill_value_kv=10.0,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "2e3b93b3", + "metadata": {}, + "source": [ + "We then add a time-dimension to our model Dataset and download knmi-data. We will have a calculation period of 3 year with daily timesteps, with a steady-state stress period with the weather of 2019 as a warmup period." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fa5fb5c3", + "metadata": {}, + "outputs": [], + "source": [ + "time = pd.date_range(\"2020\", \"2023\", freq=\"D\")\n", + "ds = nlmod.time.set_ds_time(ds, start=\"2019\", time=time, steady_start=True)\n", + "\n", + "ds.update(\n", + " nlmod.read.knmi.get_recharge(\n", + " ds, method=\"separate\", cachename=\"recharge.nc\", cachedir=ds.cachedir\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "80994109", + "metadata": {}, + "source": [ + "## Generate a simulation (sim) and groundwater flow (gwf) object\n", + "We generate a model using with all basic packages. We add drainage level at 4.0 m NAP. As the top of our model is at 6.5 m NAP this will create an unsatuated zone of about 2.5 m." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0b6b18d2", + "metadata": {}, + "outputs": [], + "source": [ + "# create simulation\n", + "sim = nlmod.sim.sim(ds)\n", + "\n", + "# create time discretisation\n", + "_ = nlmod.sim.tdis(ds, sim)\n", + "\n", + "# create groundwater flow model\n", + "gwf = nlmod.gwf.gwf(ds, sim, under_relaxation=True)\n", + "\n", + "# create ims\n", + "_ = nlmod.sim.ims(sim)\n", + "\n", + "# Create discretization\n", + "_ = nlmod.gwf.dis(ds, gwf)\n", + "\n", + "# create node property flow\n", + "_ = nlmod.gwf.npf(ds, gwf, icelltype=1)\n", + "\n", + "# creat storage\n", + "_ = nlmod.gwf.sto(ds, gwf, iconvert=1, sy=0.2, ss=1e-5)\n", + "\n", + "# Create the initial conditions package\n", + "_ = nlmod.gwf.ic(ds, gwf, starting_head=1.0)\n", + "\n", + "# Create the output control package\n", + "_ = nlmod.gwf.oc(ds, gwf)\n", + "\n", + "# set a drainage level with a resistance of 100 days in the layer that contains the drainage level\n", + "cond = ds[\"area\"] / 100.0\n", + "layer = nlmod.layers.get_layer_of_z(ds, drn_elev)\n", + "_ = nlmod.gwf.drn(ds, gwf, elev=drn_elev, cond=cond, layer=layer)" + ] + }, + { + "cell_type": "markdown", + "id": "62cbe11e", + "metadata": {}, + "source": [ + "## Run model with RCH and EVT packages\n", + "We first run the model with the Recharge (RCH) and Evapotranspiration (EVT) packages, as a reference for the model with the UZF-package." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9546d7fc", + "metadata": {}, + "outputs": [], + "source": [ + "# create recharge package\n", + "nlmod.gwf.rch(ds, gwf)\n", + "\n", + "# create evapotranspiration package\n", + "nlmod.gwf.evt(ds, gwf);" + ] + }, + { + "cell_type": "markdown", + "id": "a3f051d0", + "metadata": {}, + "source": [ + "We run this model, read the heads and get the groundwater level (the head in the highest active cell). We save the groundwater level to a variable called `gwl_rch_evt`, which we will plot later on." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5a0939fb", + "metadata": {}, + "outputs": [], + "source": [ + "# run model\n", + "nlmod.sim.write_and_run(sim, ds)\n", + "\n", + "# get heads\n", + "head_rch = nlmod.gwf.get_heads_da(ds)\n", + "\n", + "# calculate groundwater level\n", + "gwl_rch_evt = nlmod.gwf.output.get_gwl_from_wet_cells(head_rch, botm=ds[\"botm\"])" + ] + }, + { + "cell_type": "markdown", + "id": "34f1a8b6", + "metadata": {}, + "source": [ + "## Run model with UZF package\n", + "We then run the model again with the uzf-package. Before creating the uzf-package we remove the RCH- and EVT-packages.\n", + "\n", + "We choose some values for the residual water content, the saturated water content and the exponent used in the Brooks-Corey function. Other parameters are left to their defaults. The method `nlmod.gwf.uzf` will generate the UZF-package, using the variable `kv` from `ds` for the the saturated vertical hydraulic conductivity, the variable `recharge` from `ds` for the infiltration rate and the variable `evaporation` from `ds` as the potential evapotranspiration rate.\n", + "\n", + "There can be multiple layers in the unsaturated zone, just like in the saturated zone. The method `nlmod.gwf.uzf` connects the unsaturated zone cells above each other.\n", + "\n", + "Because we want to plot the water-content in the subsurface we will add some observations of the water-concent to the uzf-package. We do this by adding the optional parameter `obs_z` to `nlmod.gwf.uzf`. This will create the observations in the corresponding uzf-cells, at the requested z-values (in m NAP)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52c711d6", + "metadata": {}, + "outputs": [], + "source": [ + "gwf.remove_package(\"RCH\")\n", + "gwf.remove_package(\"EVT\")\n", + "\n", + "# create uzf package\n", + "thtr = 0.1 # residual (irreducible) water content\n", + "thts = 0.3 # saturated water content\n", + "eps = 3.5 # exponent used in the Brooks-Corey function\n", + "# add observations of the water-concent every 0.2 m, from 1 meter below the drainage-elevation to the model top\n", + "obs_z = np.arange(drn_elev-1, ds[\"top\"].max(), 0.2)[::-1]\n", + "_ = nlmod.gwf.uzf(\n", + " ds,\n", + " gwf,\n", + " thtr=thtr,\n", + " thts=thts,\n", + " thti=thtr,\n", + " eps=eps,\n", + " print_input=True,\n", + " print_flows=True,\n", + " nwavesets=100, # Modflow 6 failed and advised to increase this value from the default value of 40\n", + " obs_z=obs_z,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "b7835138", + "metadata": {}, + "source": [ + "We run the model again, and save the groundwater level to a variable called `gwl_uzf`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "39f706a5", + "metadata": {}, + "outputs": [], + "source": [ + "# run model\n", + "nlmod.sim.write_and_run(sim, ds)\n", + "\n", + "# get heads\n", + "head_rch = nlmod.gwf.get_heads_da(ds)\n", + "\n", + "# calculate groundwater level\n", + "gwl_uzf = nlmod.gwf.output.get_gwl_from_wet_cells(head_rch, botm=ds[\"botm\"])" + ] + }, + { + "cell_type": "markdown", + "id": "6938848c", + "metadata": {}, + "source": [ + "We read the water-content from the observations, and use the name of the observations to determine the layer, row, columns and z-value of each observation, which we can use in our plot." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af7128fd", + "metadata": {}, + "outputs": [], + "source": [ + "# %% \n", + "fname = os.path.join(ds.model_ws, f\"{ds.model_name}.uzf.obs.csv\")\n", + "obs = pd.read_csv(fname, index_col=0)\n", + "obs.index = pd.to_datetime(ds.time.start) + pd.to_timedelta(obs.index, \"D\")\n", + "kind, lay, row, col, z = zip(*[x.split(\"_\") for x in obs.columns])\n", + "lays = np.array([int(x) for x in lay])\n", + "rows = np.array([int(x) for x in row])\n", + "cols = np.array([int(x) for x in col])\n", + "z = np.array([float(x) for x in z])" + ] + }, + { + "cell_type": "markdown", + "id": "4d09726b", + "metadata": {}, + "source": [ + "## Compare models\n", + "We then make a plot to compare the heads in the two simulations we performed, an plot the Water Content we calculated in the UZF-calculation. We plot the water content in one vertical cell (the only) of the model. Figure layout thanks to Martin Vonk!\n", + "\n", + "The figure shows that the effect of precipitation on the groundwater level is less in summer if we also take the effect of the unsaturated zone into account (using UZF). In dry periods precipitation never reaches the groundwater level, as evaporation takes place before that." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f25f53b9", + "metadata": {}, + "outputs": [], + "source": [ + "row = 0\n", + "col = 0\n", + "\n", + "f, ax = plt.subplot_mosaic(\n", + " [[\"P\"], [\"WC\"], [\"H\"]], figsize=(10, 8), height_ratios=[2, 5, 3], sharex=True\n", + ")\n", + "\n", + "# top ax\n", + "p = ds[\"recharge\"][:, row, col].to_pandas() * 1e3\n", + "p.plot(ax=ax[\"P\"], label=\"Precipitation\", color=\"C9\")\n", + "e = ds[\"evaporation\"][:, row, col].to_pandas() * 1e3\n", + "e.plot(ax=ax[\"P\"], label=\"Evaporation\", color=\"C8\")\n", + "ax[\"P\"].set_ylabel(\"[mm/d]\")\n", + "ax[\"P\"].set_ylim(bottom=0)\n", + "ax[\"P\"].legend(loc=(0, 1), ncol=2, frameon=False)\n", + "ax[\"P\"].grid()\n", + "\n", + "# middle ax\n", + "mask = (rows == row) & (cols == col)\n", + "XY = np.meshgrid(obs.index, z[mask])\n", + "theta = obs.loc[:, mask].transpose().values\n", + "pcm = ax[\"WC\"].pcolormesh(XY[0], XY[1], theta, cmap=\"viridis_r\", vmin=thtr, vmax=thts)\n", + "ax[\"WC\"].set_ylabel(\"z [m NAP]\")\n", + "ax[\"WC\"].grid()\n", + "# set xlim, as pcolormesh increases xlim a bit\n", + "ax[\"WC\"].set_xlim(ds.time[0], ds.time[-1])\n", + "nlmod.plot.colorbar_inside(pcm, ax=ax[\"WC\"], label=r\"Moisture Content $\\theta$ [-]\")\n", + "\n", + "# bottom ax\n", + "s = gwl_rch_evt[:, row, col].to_pandas()\n", + "ax[\"H\"].plot(s.index, s.values, color=\"C1\", linestyle=\"--\", label=\"GWL RCH+EVT\")\n", + "s = gwl_uzf[:, row, col].to_pandas()\n", + "ax[\"H\"].plot(s.index, s.values, color=\"C0\", linestyle=\"-\", label=\"GWL UZF\")\n", + "ax[\"H\"].set_ylabel(\"z [m NAP]\")\n", + "ax[\"H\"].legend(loc=(0, 0.98), ncol=2, frameon=False)\n", + "ax[\"H\"].grid()\n", + "\n", + "f.tight_layout(pad=0.8)" + ] + } + ], + "metadata": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/nlmod/gwf/recharge.py b/nlmod/gwf/recharge.py index e18e5435..d25e9ac7 100644 --- a/nlmod/gwf/recharge.py +++ b/nlmod/gwf/recharge.py @@ -6,7 +6,11 @@ import xarray as xr from ..dims.grid import da_to_reclist, cols_to_reclist -from ..dims.layers import get_idomain, get_first_active_layer_from_idomain +from ..dims.layers import ( + get_idomain, + get_first_active_layer_from_idomain, + calculate_thickness, +) from ..dims.time import dataframe_to_flopy_timeseries from ..util import _get_value_from_ds_datavar @@ -55,7 +59,7 @@ def ds_to_rch(gwf, ds, mask=None, pname="rch", recharge="recharge", **kwargs): spd = da_to_reclist( ds, mask, - col1=ds[recharge].squeeze(), + col1=ds[recharge], first_active_layer=True, only_active_cells=False, ) @@ -148,7 +152,7 @@ def ds_to_evt( ds, mask, col1=surface, - col2=ds[rate].squeeze(), + col2=ds[rate], col3=depth, first_active_layer=True, only_active_cells=False, @@ -198,6 +202,8 @@ def ds_to_uzf( linear_gwet=True, unsat_etwc=False, unsat_etae=False, + obs_depth_interval=None, + obs_z=None, **kwargs, ): """Create a unsaturated zone flow package for modflow 6. This method adds uzf-cells @@ -291,6 +297,11 @@ def ds_to_uzf( If True, ET in the unsaturated zone will be simulated using a capillary pressure based formulation. Capillary pressure is calculated using the Brooks-Corey retention function. The default is False. + obs_depth_interval : float, optional + The depths at which observations of the water depth in each cell are added. The + user-specified depth must be greater than or equal to zero and less than the + thickness of GWF cellid (TOP - BOT). + The ** kwargs : dict Kwargs are passed onto flopy.mf6.ModflowGwfuzf @@ -400,11 +411,11 @@ def ds_to_uzf( if simulate_et and unsat_etae: logger.info(f"Setting root activity function (rootact) to {rootact}") - cellids = np.where(mask) + cellids_land = np.where(mask) perioddata = cols_to_reclist( ds, - cellids, + cellids_land, iuzno, finf, pet, @@ -416,6 +427,34 @@ def ds_to_uzf( cellid_column=None, ) + observations = None + # observation nodes uzf + if obs_depth_interval is not None or obs_z is not None: + cellid_per_iuzno = list(zip(*cellids)) + cellid_str = [ + str(x).replace("(", "").replace(")", "").replace(", ", "_") + for x in cellid_per_iuzno + ] + thickness = calculate_thickness(ds).data[iuzno >= 0] + obsdepths = [] + if obs_depth_interval is not None: + for i in range(nuzfcells): + depths = np.arange(obs_depth_interval, thickness[i], obs_depth_interval) + for depth in depths: + name = f"wc_{cellid_str[i]}_{depth:0.2f}" + obsdepths.append((name, "water-content", i + 1, depth)) + if obs_z is not None: + botm = ds["botm"].data[iuzno >= 0] + top = botm + thickness - landflag.data[iuzno >= 0] * surfdep / 2 + for i in range(nuzfcells): + mask = (obs_z > botm[i]) & (obs_z <= top[i]) + for z in obs_z[mask]: + depth = top[i] - z + name = f"wc_{cellid_str[i]}_{z:0.2f}" + obsdepths.append((name, "water-content", i + 1, depth)) + + observations = {ds.model_name + ".uzf.obs.csv": obsdepths} + uzf = flopy.mf6.ModflowGwfuzf( gwf, nuzfcells=nuzfcells, @@ -425,6 +464,7 @@ def ds_to_uzf( linear_gwet=linear_gwet, unsat_etwc=unsat_etwc, unsat_etae=unsat_etae, + observations=observations, **kwargs, ) diff --git a/tests/test_007_run_notebooks.py b/tests/test_007_run_notebooks.py index 985ea714..1ed78c4d 100644 --- a/tests/test_007_run_notebooks.py +++ b/tests/test_007_run_notebooks.py @@ -102,3 +102,8 @@ def test_run_notebook_15_geotop(): @pytest.mark.notebooks def test_run_notebook_16_groundwater_transport(): _run_notebook(nbdir, "16_groundwater_transport.ipynb") + + +@pytest.mark.notebooks +def test_run_notebook_17_unsaturated_zone_flow(): + _run_notebook(nbdir, "17_unsaturated_zone_flow.ipynb") diff --git a/tests/test_020_uzf.py b/tests/test_020_uzf.py index 03834201..63d0b9aa 100644 --- a/tests/test_020_uzf.py +++ b/tests/test_020_uzf.py @@ -47,4 +47,4 @@ def test_uzf_structured(): _ = nlmod.gwf.uzf(ds, gwf) # %% run - _ = nlmod.sim.write_and_run(sim, ds) + # _ = nlmod.sim.write_and_run(sim, ds) From a51eaf9abc4ac6fa32fcd887d92d8ec368ec2628 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Fri, 29 Sep 2023 17:29:22 +0200 Subject: [PATCH 03/11] Improve ds_to_rch and ds_to_evt --- nlmod/gwf/recharge.py | 33 ++++++++++++++++----------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/nlmod/gwf/recharge.py b/nlmod/gwf/recharge.py index d25e9ac7..08a903f5 100644 --- a/nlmod/gwf/recharge.py +++ b/nlmod/gwf/recharge.py @@ -43,23 +43,23 @@ def ds_to_rch(gwf, ds, mask=None, pname="rch", recharge="recharge", **kwargs): # get stress period data if ds["steady"].all(): - if "time" in ds[recharge].dims: - mask = ds[recharge].isel(time=0) != 0 - else: - mask = ds[recharge] != 0 + recharge = ds[recharge] + if "time" in recharge.dims: + recharge = recharge.isel(time=0) + mask = recharge != 0 else: rch_name_arr, rch_unique_dic = _get_unique_series(ds, recharge, pname) - recharge = "rch_name" - ds[recharge] = ds["top"].dims, rch_name_arr + ds["rch_name"] = ds["top"].dims, rch_name_arr + recharge = ds["rch_name"] if mask is not None: - mask = (ds[recharge] != "") & mask + mask = (recharge != "") & mask else: - mask = ds[recharge] != "" + mask = recharge != "" spd = da_to_reclist( ds, mask, - col1=ds[recharge], + col1=recharge, first_active_layer=True, only_active_cells=False, ) @@ -137,22 +137,21 @@ def ds_to_evt( # get stress period data if ds["steady"].all(): - if "time" in ds[rate].dims: - mask = ds[rate].isel(time=0) != 0 - else: - mask = ds[rate] != 0 + rate = ds[rate] + if "time" in rate.dims: + rate = rate.isel(time=0) + mask = rate != 0 else: evt_name_arr, evt_unique_dic = _get_unique_series(ds, rate, pname) ds["evt_name"] = ds["top"].dims, evt_name_arr - - mask = ds["evt_name"] != "" - rate = "evt_name" + rate = ds["evt_name"] + mask = rate != "" spd = da_to_reclist( ds, mask, col1=surface, - col2=ds[rate], + col2=rate, col3=depth, first_active_layer=True, only_active_cells=False, From 056b836a92254716a59e0bae5730206641e49903 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Tue, 3 Oct 2023 12:32:40 +0200 Subject: [PATCH 04/11] Add aux to rch adn evt and improve check for need for timeseries --- nlmod/dims/grid.py | 4 +-- nlmod/gwf/recharge.py | 73 +++++++++++++++++++++++++++++-------------- 2 files changed, 51 insertions(+), 26 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 049e8a4b..eff82ada 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -822,9 +822,7 @@ def lcid_to_reclist( raise ValueError("col3 is set, but col1 and/or col2 are not!") if aux is not None: - if isinstance(aux, str): - aux = [aux] - elif isinstance(aux, (int, float)): + if isinstance(aux, (str, int, float)): aux = [aux] for i_aux in aux: diff --git a/nlmod/gwf/recharge.py b/nlmod/gwf/recharge.py index 08a903f5..c83572f0 100644 --- a/nlmod/gwf/recharge.py +++ b/nlmod/gwf/recharge.py @@ -17,7 +17,9 @@ logger = logging.getLogger(__name__) -def ds_to_rch(gwf, ds, mask=None, pname="rch", recharge="recharge", **kwargs): +def ds_to_rch( + gwf, ds, mask=None, pname="rch", recharge="recharge", auxiliary=None, **kwargs +): """Convert the recharge data in the model dataset to a rch package with time series. @@ -31,6 +33,11 @@ def ds_to_rch(gwf, ds, mask=None, pname="rch", recharge="recharge", **kwargs): data array containing mask, recharge is only added where mask is True pname : str, optional package name. The default is 'rch'. + auxiliary : str or list of str + name(s) of data arrays to include as auxiliary data to reclist + recharge : str, optional + The name of the variable in ds that contains the recharge flux rate. The default + is "recharge". Returns ------- @@ -42,26 +49,27 @@ def ds_to_rch(gwf, ds, mask=None, pname="rch", recharge="recharge", **kwargs): raise ValueError("please remove nan values in recharge data array") # get stress period data - if ds["steady"].all(): + if "time" not in ds[recharge].dims or len(ds["time"]) == 1: recharge = ds[recharge] if "time" in recharge.dims: recharge = recharge.isel(time=0) - mask = recharge != 0 + mask_recharge = recharge != 0 else: rch_name_arr, rch_unique_dic = _get_unique_series(ds, recharge, pname) ds["rch_name"] = ds["top"].dims, rch_name_arr recharge = ds["rch_name"] - if mask is not None: - mask = (recharge != "") & mask - else: - mask = recharge != "" + mask_recharge = recharge != "" + + if mask is not None: + mask_recharge = mask & mask_recharge spd = da_to_reclist( ds, - mask, + mask_recharge, col1=recharge, first_active_layer=True, only_active_cells=False, + aux=auxiliary, ) # create rch package @@ -70,22 +78,35 @@ def ds_to_rch(gwf, ds, mask=None, pname="rch", recharge="recharge", **kwargs): filename=f"{gwf.name}.rch", pname=pname, fixed_cell=False, + auxiliary="CONCENTRATION" if auxiliary is not None else None, maxbound=len(spd), stress_period_data={0: spd}, **kwargs, ) + if (auxiliary is not None) and (ds.transport == 1): + logger.info("-> adding GHB to SSM sources list") + ssm_sources = list(ds.attrs["ssm_sources"]) + if rch.package_name not in ssm_sources: + ssm_sources += [rch.package_name] + ds.attrs["ssm_sources"] = ssm_sources - if ds["steady"].all(): - return rch - - # create timeseries packages - _add_time_series(rch, rch_unique_dic, ds) + if recharge.dtype == str: + # create timeseries packages + _add_time_series(rch, rch_unique_dic, ds) return rch def ds_to_evt( - gwf, ds, pname="evt", rate="evaporation", nseg=1, surface=None, depth=None, **kwargs + gwf, + ds, + mask=None, + pname="evt", + rate="evaporation", + nseg=1, + surface=None, + depth=None, + **kwargs, ): """Convert the evaporation data in the model dataset to a evt package with time series. @@ -96,8 +117,13 @@ def ds_to_evt( groundwater flow model. ds : xr.DataSet dataset containing relevant model grid information + mask : xr.DataArray + data array containing mask, evt is only added where mask is True pname : str, optional package name. The default is 'evt'. + rate : str, optional + The name of the variable in ds that contains the maximum ET flux rate. The + default is "evaporation". nseg : int, optional number of ET segments. Only 1 is supported for now. The default is 1. surface : str, float or xr.DataArray, optional @@ -136,20 +162,23 @@ def ds_to_evt( raise ValueError("please remove nan values in evaporation data array") # get stress period data - if ds["steady"].all(): + if "time" not in ds[rate].dims or len(ds["time"]) == 1: rate = ds[rate] if "time" in rate.dims: rate = rate.isel(time=0) - mask = rate != 0 + mask_rate = rate != 0 else: evt_name_arr, evt_unique_dic = _get_unique_series(ds, rate, pname) ds["evt_name"] = ds["top"].dims, evt_name_arr rate = ds["evt_name"] - mask = rate != "" + mask_rate = rate != "" + + if mask is not None: + mask_rate = mask & mask_rate spd = da_to_reclist( ds, - mask, + mask_rate, col1=surface, col2=rate, col3=depth, @@ -169,11 +198,9 @@ def ds_to_evt( **kwargs, ) - if ds["steady"].all(): - return evt - - # create timeseries packages - _add_time_series(evt, evt_unique_dic, ds) + if rate.dtype == str: + # create timeseries packages + _add_time_series(evt, evt_unique_dic, ds) return evt From d22fdc7ab80154a782e8c505ec52d36d617b914d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Tue, 3 Oct 2023 13:55:55 +0200 Subject: [PATCH 05/11] solve bug where timeseries were not written --- nlmod/gwf/recharge.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/nlmod/gwf/recharge.py b/nlmod/gwf/recharge.py index c83572f0..f4b56097 100644 --- a/nlmod/gwf/recharge.py +++ b/nlmod/gwf/recharge.py @@ -49,7 +49,8 @@ def ds_to_rch( raise ValueError("please remove nan values in recharge data array") # get stress period data - if "time" not in ds[recharge].dims or len(ds["time"]) == 1: + use_ts = "time" in ds[recharge].dims and len(ds["time"]) > 1 + if not use_ts: recharge = ds[recharge] if "time" in recharge.dims: recharge = recharge.isel(time=0) @@ -90,7 +91,7 @@ def ds_to_rch( ssm_sources += [rch.package_name] ds.attrs["ssm_sources"] = ssm_sources - if recharge.dtype == str: + if use_ts: # create timeseries packages _add_time_series(rch, rch_unique_dic, ds) @@ -162,7 +163,8 @@ def ds_to_evt( raise ValueError("please remove nan values in evaporation data array") # get stress period data - if "time" not in ds[rate].dims or len(ds["time"]) == 1: + use_ts = "time" in ds[rate].dims and len(ds["time"]) > 1 + if not use_ts: rate = ds[rate] if "time" in rate.dims: rate = rate.isel(time=0) @@ -198,7 +200,7 @@ def ds_to_evt( **kwargs, ) - if rate.dtype == str: + if use_ts: # create timeseries packages _add_time_series(evt, evt_unique_dic, ds) From f954eb104ab6bd7f1be3d24df091cd400fbf8af7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Wed, 4 Oct 2023 21:36:51 +0200 Subject: [PATCH 06/11] Improve wel_from_df and maw_from_df Distribute k based on kh-value in wel_from_df Add multiplier to wel_from_df, so q can be a timeseries Add aux to maw_from_df Improve speed by using gdf_to_grid instead of slow gwf.modelgrid.intersect becuase it is fast, remove tqdm maw-cells can be added to top layer as well (was not possible before) --- nlmod/gwf/wells.py | 213 +++++++++++++++++++++++++++++++-------------- 1 file changed, 148 insertions(+), 65 deletions(-) diff --git a/nlmod/gwf/wells.py b/nlmod/gwf/wells.py index 5ba7c9a8..65037154 100644 --- a/nlmod/gwf/wells.py +++ b/nlmod/gwf/wells.py @@ -2,7 +2,10 @@ import flopy as fp import numpy as np -from tqdm import tqdm +import pandas as pd +import geopandas as gpd + +from ..dims.grid import gdf_to_grid logger = logging.getLogger(__name__) @@ -17,46 +20,44 @@ def wel_from_df( Q="Q", aux=None, boundnames=None, + ds=None, + auxmultname="multiplier", **kwargs, ): + if aux is None: + aux = [] + if not isinstance(aux, list): + aux = [aux] + + df = _add_cellid(df, ds=ds, gwf=gwf, x=x, y=y) + multipliers = _get_layer_multiplier_for_wells(df, top, botm, ds=ds, gwf=gwf) + # collect data well_lrcd = [] - - for wnam, irow in tqdm(df.iterrows(), total=df.index.size, desc="Adding WELs"): - try: - cid1 = gwf.modelgrid.intersect(irow[x], irow[y], irow[top], forgive=False) - cid2 = gwf.modelgrid.intersect(irow[x], irow[y], irow[botm], forgive=False) - except Exception: - logger.warning( - f"Warning! well {wnam} outside of model domain! ({irow[x]}, {irow[y]})" - ) - continue - kb = cid2[0] - if len(cid1) == 2: - kt, icell2d = cid1 - idomain_mask = gwf.modelgrid.idomain[kt : kb + 1, icell2d] > 0 - elif len(cid1) == 3: - kt, i, j = cid1 - idomain_mask = gwf.modelgrid.idomain[kt : kb + 1, i, j] > 0 - # mask only active cells - wlayers = np.arange(kt, kb + 1)[idomain_mask] - if len(wlayers) > 0: - for k in wlayers: - if len(cid1) == 2: - wdata = [(k, icell2d), irow[Q] / len(wlayers)] - elif len(cid1) == 3: - wdata = [(k, i, j), irow[Q] / len(wlayers)] - - if aux is not None: - if not isinstance(aux, list): - aux = [aux] - for iaux in aux: - wdata.append(irow[iaux]) - if boundnames is not None: - wdata.append(irow[boundnames]) - well_lrcd.append(wdata) - else: - logger.warning(f"Well '{wnam}' not added, no active screened layers.") + for index, irow in df.iterrows(): + wlayers = np.where(multipliers[index] > 0)[0] + for k in wlayers: + multiplier = multipliers[index][k] + q = irow[Q] + if auxmultname is None: + q = q * multiplier + if isinstance(irow["cellid"], int): + # vertex grid + cellid = (k, irow["cellid"]) + else: + # structured grid + cellid = (k, irow["cellid"][0], irow["cellid"][1]) + wdata = [cellid, q] + for iaux in aux: + wdata.append(irow[iaux]) + if auxmultname is not None: + wdata.append(multiplier) + if boundnames is not None: + wdata.append(irow[boundnames]) + well_lrcd.append(wdata) + + if auxmultname is not None: + aux.append(auxmultname) wel_spd = {0: well_lrcd} @@ -65,6 +66,7 @@ def wel_from_df( stress_period_data=wel_spd, auxiliary=aux, boundnames=boundnames is not None, + auxmultname=auxmultname, **kwargs, ) @@ -82,16 +84,25 @@ def maw_from_df( rw="rw", condeqn="THIEM", strt=0.0, + aux=None, boundnames=None, + ds=None, **kwargs, ): + if aux is None: + aux = [] + if not isinstance(aux, list): + aux = [aux] + + df = _add_cellid(df, ds=ds, gwf=gwf, x=x, y=y) + multipliers = _get_layer_multiplier_for_wells(df, top, botm, ds=ds, gwf=gwf) + maw_pakdata = [] maw_conndata = [] maw_perdata = [] iw = 0 - for index in tqdm(df.index, total=df.index.size, desc="Adding MAWs"): - irow = df.loc[index] + for index, irow in df.iterrows(): try: cid1 = gwf.modelgrid.intersect(irow[x], irow[y], irow[top], forgive=False) cid2 = gwf.modelgrid.intersect(irow[x], irow[y], irow[botm], forgive=False) @@ -110,8 +121,11 @@ def maw_from_df( wlayers = np.arange(kt, kb + 1)[idomain_mask] + wlayers = np.where(multipliers[index] > 0)[0] # pakdata = [iw, irow[rw], irow[botm], strt, condeqn, len(wlayers)] + for iaux in aux: + pakdata.append(irow[iaux]) if boundnames is not None: pakdata.append(irow[boundnames]) maw_pakdata.append(pakdata) @@ -119,38 +133,41 @@ def maw_from_df( maw_perdata.append([iw, "RATE", irow[Q]]) for iwellpart, k in enumerate(wlayers): - if gwf.modelgrid.grid_type == "vertex": - laytop = gwf.modelgrid.botm[k - 1, icell2d] - laybot = gwf.modelgrid.botm[k, icell2d] - # - mawdata = [ - iw, - iwellpart, - (k, icell2d), - np.min([irow[top], laytop]), - laybot, - 0.0, - 0.0, - ] - elif gwf.modelgrid.grid_type == "structured": - laytop = gwf.modelgrid.botm[k - 1, i, j] - laybot = gwf.modelgrid.botm[k, i, j] - # - mawdata = [ - iw, - iwellpart, - (k, i, j), - np.min([irow[top], laytop]), - laybot, - 0.0, - 0.0, - ] + if k == 0: + laytop = gwf.modelgrid.top + else: + laytop = gwf.modelgrid.botm[k - 1] + laybot = gwf.modelgrid.botm[k] + + if isinstance(irow["cellid"], int): + # vertex grid + cellid = (k, irow["cellid"]) + laytop = laytop[irow["cellid"]] + laybot = laybot[irow["cellid"]] + else: + # structured grid + cellid = (k, irow["cellid"][0], irow["cellid"][1]) + laytop = laytop[irow["cellid"][0], irow["cellid"][1]] + laybot = laybot[irow["cellid"][0], irow["cellid"][1]] + scrn_top = np.min([irow[top], laytop]) + scrn_bot = np.max([irow[botm], laybot]) + + mawdata = [ + iw, + iwellpart, + cellid, + scrn_top, + scrn_bot, + 0.0, + 0.0, + ] maw_conndata.append(mawdata) iw += 1 maw = fp.mf6.ModflowGwfmaw( gwf, nmawwells=iw, + auxiliary=aux, boundnames=boundnames is not None, packagedata=maw_pakdata, connectiondata=maw_conndata, @@ -159,3 +176,69 @@ def maw_from_df( ) return maw + + +def _add_cellid(df, ds=None, gwf=None, x=None, y=None): + if not isinstance(df, gpd.GeoDataFrame): + df = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df[x], df[y])) + if "cellid" not in df.columns: + df = gdf_to_grid(df, gwf if ds is None else ds) + return df + + +def _get_layer_multiplier_for_wells(df, top, botm, ds=None, gwf=None): + # get required data either from gwf or ds + if ds is None: + ml_top = gwf.dis.top.array + ml_bot = gwf.dis.botm.array + kh = gwf.npf.k.array + layer = range(gwf.dis.nlay.array) + else: + ml_top = ds["top"].data + ml_bot = ds["botm"].data + kh = ds["kh"].data + layer = ds.layer + + multipliers = dict() + for index, irow in df.iterrows(): + multipliers[index] = _get_layer_multiplier_for_well( + irow["cellid"], irow[top], irow[botm], ml_top, ml_bot, kh + ) + + if (multipliers[index] == 0).all(): + logger.warning(f"No layers found for well {index}") + multipliers = pd.DataFrame(multipliers, index=layer, columns=df.index) + return multipliers + + +def _get_layer_multiplier_for_well(cid, well_top, well_bot, ml_top, ml_bot, ml_kh): + """Get a factor for each layer that a well""" + # keep the tops and botms of the cell where the well is in + ml_top_cid = ml_top[cid].copy() + if isinstance(cid, int): + ml_bot_cid = ml_bot[:, cid].copy() + ml_kh_cid = ml_kh[:, cid].copy() + else: + ml_bot_cid = ml_bot[:, cid[0], cid[1]].copy() + ml_kh_cid = ml_kh[:, cid[0], cid[1]].copy() + ml_top_cid = np.array([ml_top_cid] + list(ml_bot_cid[:-1])) + + # only keep the part of layers along the well filter + ml_top_cid[ml_top_cid > well_top] = well_top + ml_top_cid[ml_top_cid < well_bot] = well_bot + ml_bot_cid[ml_bot_cid > well_top] = well_top + ml_bot_cid[ml_bot_cid < well_bot] = well_bot + + # calculate remaining kd along the well filter + kd = ml_kh_cid * (ml_top_cid - ml_bot_cid) + mask = kd < 0 + if np.any(mask): + logger.warning("There are negative thicknesses at cellid {cid}") + kd[mask] = 0 + if (kd == 0).all(): + # the well does not cross any of the layers. Just return an array of zeros. + multiplier = kd + else: + # devide by the total kd to get a factor + multiplier = kd / kd.sum() + return multiplier From 802b5eea44814ba612d9bdd5bad61bac208e6bd1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Wed, 4 Oct 2023 22:04:11 +0200 Subject: [PATCH 07/11] Remove some metadata from notebooks --- docs/examples/00_model_from_scratch.ipynb | 8 +------- docs/examples/02_surface_water.ipynb | 8 +------- docs/examples/09_schoonhoven.ipynb | 10 ++-------- docs/examples/11_grid_rotation.ipynb | 8 +------- docs/examples/14_stromingen_example.ipynb | 8 +------- 5 files changed, 6 insertions(+), 36 deletions(-) diff --git a/docs/examples/00_model_from_scratch.ipynb b/docs/examples/00_model_from_scratch.ipynb index b043358b..40579e4f 100644 --- a/docs/examples/00_model_from_scratch.ipynb +++ b/docs/examples/00_model_from_scratch.ipynb @@ -321,14 +321,8 @@ } ], "metadata": { - "kernelspec": { - "display_name": "artesia", - "language": "python", - "name": "python3" - }, "language_info": { - "name": "python", - "version": "3.10.12" + "name": "python" } }, "nbformat": 4, diff --git a/docs/examples/02_surface_water.ipynb b/docs/examples/02_surface_water.ipynb index 8e705bd7..22a1217e 100644 --- a/docs/examples/02_surface_water.ipynb +++ b/docs/examples/02_surface_water.ipynb @@ -798,14 +798,8 @@ } ], "metadata": { - "kernelspec": { - "display_name": "artesia", - "language": "python", - "name": "python3" - }, "language_info": { - "name": "python", - "version": "3.10.12" + "name": "python" } }, "nbformat": 4, diff --git a/docs/examples/09_schoonhoven.ipynb b/docs/examples/09_schoonhoven.ipynb index 530f0852..266d28b2 100644 --- a/docs/examples/09_schoonhoven.ipynb +++ b/docs/examples/09_schoonhoven.ipynb @@ -338,7 +338,7 @@ " \"identificatie\"\n", "].isin(ids_oude_haven)\n", "lakes = bgt_grid[mask].copy()\n", - "lakes[\"name\"] = \"\"\n", + "lakes[\"name\"] = \"\"\n", "lakes.loc[lakes[\"identificatie\"].isin(ids_grote_gracht), \"name\"] = \"grote gracht\"\n", "lakes.loc[lakes[\"identificatie\"].isin(ids_oude_haven), \"name\"] = \"oude haven\"\n", "bgt_grid = bgt_grid[~mask]\n" @@ -827,14 +827,8 @@ } ], "metadata": { - "kernelspec": { - "display_name": "nlmod", - "language": "python", - "name": "python3" - }, "language_info": { - "name": "python", - "version": "3.9.4" + "name": "python" } }, "nbformat": 4, diff --git a/docs/examples/11_grid_rotation.ipynb b/docs/examples/11_grid_rotation.ipynb index 852f2efb..e522cc55 100644 --- a/docs/examples/11_grid_rotation.ipynb +++ b/docs/examples/11_grid_rotation.ipynb @@ -358,14 +358,8 @@ } ], "metadata": { - "kernelspec": { - "display_name": "artesia", - "language": "python", - "name": "python3" - }, "language_info": { - "name": "python", - "version": "3.10.12" + "name": "python" } }, "nbformat": 4, diff --git a/docs/examples/14_stromingen_example.ipynb b/docs/examples/14_stromingen_example.ipynb index 41f4ab96..0dea2056 100644 --- a/docs/examples/14_stromingen_example.ipynb +++ b/docs/examples/14_stromingen_example.ipynb @@ -330,14 +330,8 @@ } ], "metadata": { - "kernelspec": { - "display_name": "artesia", - "language": "python", - "name": "python3" - }, "language_info": { - "name": "python", - "version": "3.10.12" + "name": "python" } }, "nbformat": 4, From 0cdbbdb62e927f03d79e872b7fdfc68be33c8985 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Wed, 4 Oct 2023 22:31:31 +0200 Subject: [PATCH 08/11] Fix bugs found by tests, and add extra tests --- nlmod/gwf/wells.py | 5 +++++ tests/test_010_wells.py | 24 ++++++++++++++++++++++++ 2 files changed, 29 insertions(+) diff --git a/nlmod/gwf/wells.py b/nlmod/gwf/wells.py index 65037154..5905cb99 100644 --- a/nlmod/gwf/wells.py +++ b/nlmod/gwf/wells.py @@ -61,6 +61,9 @@ def wel_from_df( wel_spd = {0: well_lrcd} + if len(aux) == 0: + aux = None + wel = fp.mf6.ModflowGwfwel( gwf, stress_period_data=wel_spd, @@ -164,6 +167,8 @@ def maw_from_df( maw_conndata.append(mawdata) iw += 1 + if len(aux) == 0: + aux = None maw = fp.mf6.ModflowGwfmaw( gwf, nmawwells=iw, diff --git a/tests/test_010_wells.py b/tests/test_010_wells.py index 26a2adc9..bd07b143 100644 --- a/tests/test_010_wells.py +++ b/tests/test_010_wells.py @@ -50,6 +50,15 @@ def test_wel_from_df(): nlmod.gwf.wells.wel_from_df(wells, gwf) +def test_wel_from_df_no_multiplier(): + wells = pd.DataFrame(columns=["x", "y", "top", "botm", "Q"], index=range(2)) + wells.loc[0] = 100, -50, -5, -10, -100.0 + wells.loc[1] = 200, 150, -20, -30, -300.0 + + sim, gwf = get_sim_and_gwf() + nlmod.gwf.wells.wel_from_df(wells, gwf, auxmultname=None) + + def test_maw_from_df(): wells = pd.DataFrame(columns=["x", "y", "top", "botm", "rw", "Q"], index=range(2)) wells.loc[0] = 100, -50, -5, -10, 0.1, -100.0 @@ -59,6 +68,21 @@ def test_maw_from_df(): nlmod.gwf.wells.maw_from_df(wells, gwf) +def test_wel_from_df_transient(): + time = pd.date_range("2000", "2001", freq="MS") + ds = get_model_ds(start="1990", time=time) + + Q = pd.DataFrame(-100.0, index=time, columns=["Q1", "Q2"]) + wells = pd.DataFrame(columns=["x", "y", "top", "botm", "Q"], index=range(2)) + wells.loc[0] = 100, -50, -5, -10, "Q1" + wells.loc[1] = 200, 150, -20, -30, "Q2" + + sim, gwf = get_sim_and_gwf() + wel = nlmod.gwf.wells.wel_from_df(wells, gwf) + + nlmod.time.dataframe_to_flopy_timeseries(Q, ds, package=wel) + + def test_maw_from_df_transient(): time = pd.date_range("2000", "2001", freq="MS") ds = get_model_ds(start="1990", time=time) From d6041c14d3178fda8a2320a8d4e2f31a2e959c4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Wed, 4 Oct 2023 23:06:23 +0200 Subject: [PATCH 09/11] Fix codacy issues and some spelling in the uzf-notebook --- docs/examples/17_unsaturated_zone_flow.ipynb | 12 ++++++------ nlmod/dims/grid.py | 4 +--- nlmod/gwf/recharge.py | 4 +++- nlmod/gwf/wells.py | 2 +- 4 files changed, 11 insertions(+), 11 deletions(-) diff --git a/docs/examples/17_unsaturated_zone_flow.ipynb b/docs/examples/17_unsaturated_zone_flow.ipynb index 57aa3c80..31a19a77 100644 --- a/docs/examples/17_unsaturated_zone_flow.ipynb +++ b/docs/examples/17_unsaturated_zone_flow.ipynb @@ -109,7 +109,7 @@ "id": "cdb20bc9", "metadata": {}, "source": [ - "As you can see, there are some NaN-values in the hydaulic conductivities (`kh` and `kv`). These will be filled when making a model Dataset, usinge fill-values and anisotropy values. See the info-messages after the commands below." + "As you can see, there are some NaN-values in the hydaulic conductivities (`kh` and `kv`). These will be filled when making a model Dataset, using fill-values and anisotropy values. See the info-messages after the commands below." ] }, { @@ -258,11 +258,11 @@ "## Run model with UZF package\n", "We then run the model again with the uzf-package. Before creating the uzf-package we remove the RCH- and EVT-packages.\n", "\n", - "We choose some values for the residual water content, the saturated water content and the exponent used in the Brooks-Corey function. Other parameters are left to their defaults. The method `nlmod.gwf.uzf` will generate the UZF-package, using the variable `kv` from `ds` for the the saturated vertical hydraulic conductivity, the variable `recharge` from `ds` for the infiltration rate and the variable `evaporation` from `ds` as the potential evapotranspiration rate.\n", + "We choose some values for the residual water content, the saturated water content and the exponent used in the Brooks-Corey function. Other parameters are left to their defaults. The method `nlmod.gwf.uzf` will generate the UZF-package, using the variable `kv` from `ds` for the saturated vertical hydraulic conductivity, the variable `recharge` from `ds` for the infiltration rate and the variable `evaporation` from `ds` as the potential evapotranspiration rate.\n", "\n", "There can be multiple layers in the unsaturated zone, just like in the saturated zone. The method `nlmod.gwf.uzf` connects the unsaturated zone cells above each other.\n", "\n", - "Because we want to plot the water-content in the subsurface we will add some observations of the water-concent to the uzf-package. We do this by adding the optional parameter `obs_z` to `nlmod.gwf.uzf`. This will create the observations in the corresponding uzf-cells, at the requested z-values (in m NAP)." + "Because we want to plot the water content in the subsurface we will add some observations of the water concent to the uzf-package. We do this by adding the optional parameter `obs_z` to `nlmod.gwf.uzf`. This will create the observations in the corresponding uzf-cells, at the requested z-values (in m NAP)." ] }, { @@ -279,7 +279,7 @@ "thtr = 0.1 # residual (irreducible) water content\n", "thts = 0.3 # saturated water content\n", "eps = 3.5 # exponent used in the Brooks-Corey function\n", - "# add observations of the water-concent every 0.2 m, from 1 meter below the drainage-elevation to the model top\n", + "# add observations of the water concent every 0.2 m, from 1 meter below the drainage-elevation to the model top\n", "obs_z = np.arange(drn_elev-1, ds[\"top\"].max(), 0.2)[::-1]\n", "_ = nlmod.gwf.uzf(\n", " ds,\n", @@ -325,7 +325,7 @@ "id": "6938848c", "metadata": {}, "source": [ - "We read the water-content from the observations, and use the name of the observations to determine the layer, row, columns and z-value of each observation, which we can use in our plot." + "We read the water content from the observations, and use the name of the observations to determine the layer, row, columns and z-value of each observation, which we can use in our plot." ] }, { @@ -352,7 +352,7 @@ "metadata": {}, "source": [ "## Compare models\n", - "We then make a plot to compare the heads in the two simulations we performed, an plot the Water Content we calculated in the UZF-calculation. We plot the water content in one vertical cell (the only) of the model. Figure layout thanks to Martin Vonk!\n", + "We then make a plot to compare the heads in the two simulations we performed, and plot the water content we calculated in the UZF-calculation, and added observations for. We plot the water content in one vertical cell (the only) of the model. Figure layout thanks to Martin Vonk!\n", "\n", "The figure shows that the effect of precipitation on the groundwater level is less in summer if we also take the effect of the unsaturated zone into account (using UZF). In dry periods precipitation never reaches the groundwater level, as evaporation takes place before that." ] diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index eff82ada..3402057e 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -727,9 +727,7 @@ def lrc_to_reclist( raise ValueError("col3 is set, but col1 and/or col2 are not!") if aux is not None: - if isinstance(aux, str): - aux = [aux] - elif isinstance(aux, (int, float)): + if isinstance(aux, (str, int, float)): aux = [aux] for i_aux in aux: diff --git a/nlmod/gwf/recharge.py b/nlmod/gwf/recharge.py index f4b56097..4bf9e659 100644 --- a/nlmod/gwf/recharge.py +++ b/nlmod/gwf/recharge.py @@ -485,6 +485,8 @@ def ds_to_uzf( uzf = flopy.mf6.ModflowGwfuzf( gwf, + filename=f"{gwf.name}.uzf", + pname=pname, nuzfcells=nuzfcells, packagedata=packagedata, perioddata={0: perioddata}, @@ -579,7 +581,7 @@ def _add_time_series(package, rch_unique_dic, ds): df = pd.DataFrame(rch_unique_dic, index=ds.time) if df.isna().any(axis=None): # make sure there are no NaN's, as otherwise they will be filled by zeros later - raise (Exception("There cannot be nan's in the dictionary")) + raise (ValueError("There cannot be nan's in the DataFrame")) # set the first value for the start-time as well df.loc[pd.to_datetime(ds.time.start)] = df.iloc[0] # combine the values with the start of each period, and ste the last value to 0.0 diff --git a/nlmod/gwf/wells.py b/nlmod/gwf/wells.py index 5905cb99..9c757273 100644 --- a/nlmod/gwf/wells.py +++ b/nlmod/gwf/wells.py @@ -204,7 +204,7 @@ def _get_layer_multiplier_for_wells(df, top, botm, ds=None, gwf=None): kh = ds["kh"].data layer = ds.layer - multipliers = dict() + multipliers = {} for index, irow in df.iterrows(): multipliers[index] = _get_layer_multiplier_for_well( irow["cellid"], irow[top], irow[botm], ml_top, ml_bot, kh From 330ec62b3665fdf45be4f1e322e6774f515cefc5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Thu, 5 Oct 2023 20:05:53 +0200 Subject: [PATCH 10/11] Process comments and add docstrings to wells.py --- nlmod/gwf/gwf.py | 8 +- nlmod/gwf/recharge.py | 44 ++++++---- nlmod/gwf/wells.py | 185 ++++++++++++++++++++++++++++++++++++++++-- 3 files changed, 206 insertions(+), 31 deletions(-) diff --git a/nlmod/gwf/gwf.py b/nlmod/gwf/gwf.py index 8653fbe0..2ce266db 100644 --- a/nlmod/gwf/gwf.py +++ b/nlmod/gwf/gwf.py @@ -395,7 +395,6 @@ def ghb( ghb = flopy.mf6.ModflowGwfghb( gwf, auxiliary="CONCENTRATION" if auxiliary is not None else None, - print_input=True, maxbound=len(ghb_rec), stress_period_data=ghb_rec, save_flows=True, @@ -472,7 +471,6 @@ def drn( if len(drn_rec) > 0: drn = flopy.mf6.ModflowGwfdrn( gwf, - print_input=True, maxbound=len(drn_rec), stress_period_data=drn_rec, save_flows=True, @@ -552,7 +550,6 @@ def riv( if len(riv_rec) > 0: riv = flopy.mf6.ModflowGwfriv( gwf, - print_input=True, maxbound=len(riv_rec), stress_period_data=riv_rec, auxiliary="CONCENTRATION" if auxiliary is not None else None, @@ -795,7 +792,6 @@ def surface_drain_from_ds(ds, gwf, resistance, elev="ahn", pname="drn", **kwargs drn = flopy.mf6.ModflowGwfdrn( gwf, pname=pname, - print_input=True, maxbound=len(drn_rec), stress_period_data={0: drn_rec}, save_flows=True, @@ -874,9 +870,9 @@ def uzf(ds, gwf, pname="uzf", **kwargs): logger.info("creating mf6 UZF") # create uzf package - evt = recharge.ds_to_uzf(gwf, ds, pname=pname, **kwargs) + uzf = recharge.ds_to_uzf(gwf, ds, pname=pname, **kwargs) - return evt + return uzf def _set_record(out, budget, output="head"): diff --git a/nlmod/gwf/recharge.py b/nlmod/gwf/recharge.py index 4bf9e659..07e35cfc 100644 --- a/nlmod/gwf/recharge.py +++ b/nlmod/gwf/recharge.py @@ -280,8 +280,9 @@ def ds_to_uzf( landflag : xr.DataArray, optional A DataArray with integer values, set to one for land surface cells indicating that boundary conditions can be applied and data can be specified in the PERIOD - block. A value of 0 specifies a non-land surface cell. Landflag is determined - from ds when it is None. The default is None. + block. A value of 0 specifies a non-land surface cell. Landflag is set to one + for the most upper active layer in each vertical column (determined form ds) + when it is None. The default is None. finf :float, str or array-like The applied infiltration rate of the UZF cell (:math:`LT^{-1}`). When passed as string, the array is obtained from ds. The default is "recharge". @@ -289,27 +290,28 @@ def ds_to_uzf( The potential evapotranspiration rate of the UZF cell and specified GWF cell. Evapotranspiration is first removed from the unsaturated zone and any remaining potential evapotranspiration is applied to the saturated zone. only - used if SIMULATE_ET. When passed as string, the array is obtained from ds. The - default is "evaporation". + used if simulate_et=True. When passed as string, the array is obtained from ds. + The default is "evaporation". extdp : float, optional Value that defines the evapotranspiration extinction depth of the UZF cell, in meters below the top of the model. Set to 2.0 meter when None. The default is None. extwc : float, optional The evapotranspiration extinction water content of the UZF cell. Only used if - SIMULATE_ET and UNSAT_ETWC. Set to thtr when None. The default is None. + simulate_et=True and unsat_etwc=True. Set to thtr when None. The default is + None. ha : float, optional - The air entry potential (head) of the UZF cell. Only used if SIMULATE_ET and - UNSAT_ETAE. Set to 0.0 when None. The default is None. + The air entry potential (head) of the UZF cell. Only used if simulate_et=True + and unsat_etae=True. Set to 0.0 when None. The default is None. hroot : float, optional - The root potential (head) of the UZF cell. Only used if SIMULATE_ET and - UNSAT_ETAE. Set to 0.0 when None. The default is None. + The root potential (head) of the UZF cell. Only used if simulate_et=True and + unsat_etae=True. Set to 0.0 when None. The default is None. rootact : float, optional the root activity function of the UZF cell. ROOTACT is the length of roots in a given volume of soil divided by that volume. Values range from 0 to about 3 :math:`cm^{-2}`, depending on the plant community and its stage of development. - Only used if SIMULATE_ET and UNSAT_ETAE. Set to 0.0 when None. The default is - None. + Only used if simulate_et=True and unsat_etae=True. Set to 0.0 when None. The + default is None. simulate_et : bool, optional If True, ET in the unsaturated (UZF) and saturated zones (GWF) will be simulated. The default is True. @@ -326,10 +328,13 @@ def ds_to_uzf( based formulation. Capillary pressure is calculated using the Brooks-Corey retention function. The default is False. obs_depth_interval : float, optional - The depths at which observations of the water depth in each cell are added. The - user-specified depth must be greater than or equal to zero and less than the - thickness of GWF cellid (TOP - BOT). - The + The depths at which observations of the water content in each cell are added. + When not None, this creates a CSV output file with water content at different + z-coordinates in each UZF cell. The default is None. + obs_z : array-like, optional + The z-coordinate at which observations of the water content in each cell are + added. When not None, this creates a CSV output file with water content at fixes + z-coordinates in each UZF cell. The default is None. ** kwargs : dict Kwargs are passed onto flopy.mf6.ModflowGwfuzf @@ -464,6 +469,10 @@ def ds_to_uzf( for x in cellid_per_iuzno ] thickness = calculate_thickness(ds).data[iuzno >= 0] + # account for surfdep, as this decreases the height of the top of the upper cell + # otherwise modflow may return an error + thickness = thickness - landflag.data[iuzno >= 0] * surfdep / 2 + obsdepths = [] if obs_depth_interval is not None: for i in range(nuzfcells): @@ -473,7 +482,8 @@ def ds_to_uzf( obsdepths.append((name, "water-content", i + 1, depth)) if obs_z is not None: botm = ds["botm"].data[iuzno >= 0] - top = botm + thickness - landflag.data[iuzno >= 0] * surfdep / 2 + + top = botm + thickness for i in range(nuzfcells): mask = (obs_z > botm[i]) & (obs_z <= top[i]) for z in obs_z[mask]: @@ -584,7 +594,7 @@ def _add_time_series(package, rch_unique_dic, ds): raise (ValueError("There cannot be nan's in the DataFrame")) # set the first value for the start-time as well df.loc[pd.to_datetime(ds.time.start)] = df.iloc[0] - # combine the values with the start of each period, and ste the last value to 0.0 + # combine the values with the start of each period, and set the last value to 0.0 df = df.sort_index().shift(-1).fillna(0.0) dataframe_to_flopy_timeseries( diff --git a/nlmod/gwf/wells.py b/nlmod/gwf/wells.py index 9c757273..d9db0b81 100644 --- a/nlmod/gwf/wells.py +++ b/nlmod/gwf/wells.py @@ -24,6 +24,52 @@ def wel_from_df( auxmultname="multiplier", **kwargs, ): + """ + Add a Well (WEL) package based on input from a (Geo)DataFrame. + + Parameters + ---------- + df : pd.DataFrame or gpd.GeoDataFrame + A (Geo)DataFrame containing the properties of the wells. + gwf : flopy ModflowGwf + Groundwaterflow object to add the wel-package to. + x : str, optional + The column in df that contains the x-coordinate of the well. Only used when df + is a DataFrame. The default is 'x'. + y : str, optional + The column in df that contains the y-coordinate of the well. Only used when df + is a DataFrame. The default is 'y'. + top : str + The column in df that contains the z-coordinate of the top of the well screen. + The defaults is 'top'. + botm : str + The column in df that contains the z-coordinate of the bottom of the well + screen. The defaults is 'botm'. + Q : str, optional + The column in df containing the well discharge. This column can contain floats, + or strings belonging to timeseries added later. The default is "Q". + aux : str of list of str, optional + The column(s) in df that contain auxiliary variables. The default is None. + boundnames : str, optional + THe column in df thet . The default is None. + ds : xarray.Dataset + Dataset with model data. Needed to determine cellid when grid-rotation is used. + The default is None. + auxmultname : str, optional + The name of the auxiliary varibale that contains the multiplication factors to + distribute the well discharge over different layers. When auxmultname is None, + this auxiliary variable will not be added, and Q is multiplied by these factors + directly. auxmultname cannot be None when df[Q] contains strings (the names of + timeseries). The default is "multiplier". + **kwargs : dict + kwargs are passed to flopy.mf6.ModflowGwfwel. + + Returns + ------- + wel : flopy.mf6.ModflowGwfwel + wel package. + + """ if aux is None: aux = [] if not isinstance(aux, list): @@ -92,6 +138,54 @@ def maw_from_df( ds=None, **kwargs, ): + """ + Add a Multi-AquiferWell (MAW) package based on input from a (Geo)DataFrame. + + Parameters + ---------- + df : pd.DataFrame or gpd.GeoDataFrame + A (Geo)DataFrame containing the properties of the wells. + gwf : flopy ModflowGwf + Groundwaterflow object to add the wel-package to. + x : str, optional + The column in df that contains the x-coordinate of the well. Only used when df + is a DataFrame. The default is 'x'. + y : str, optional + The column in df that contains the y-coordinate of the well. Only used when df + is a DataFrame. The default is 'y'. + top : str + The column in df that contains the z-coordinate of the top of the well screen. + The defaults is 'top'. + botm : str + The column in df that contains the z-coordinate of the bottom of the well + screen. The defaults is 'botm'. + Q : str, optional + The column in df containing the well discharge. This column can contain floats, + or strings belonging to timeseries added later. The default is "Q". + rw : str, optional + The column in df that contains the radius for the multi-aquifer well. The + default is "rw". + condeqn : str, optional + String that defines the conductance equation that is used to calculate the + saturated conductance for the multi-aquifer well. The default is "THIEM". + strt : float, optional + The starting head for the multi-aquifer well. The default is 0.0. + aux : str of list of str, optional + The column(s) in df that contain auxiliary variables. The default is None. + boundnames : str, optional + THe column in df thet . The default is None. + ds : xarray.Dataset + Dataset with model data. Needed to determine cellid when grid-rotation is used. + The default is None. + **kwargs : TYPE + Kwargs are passed onto ModflowGwfmaw. + + Returns + ------- + wel : flopy.mf6.ModflowGwfmaw + maw package. + + """ if aux is None: aux = [] if not isinstance(aux, list): @@ -183,7 +277,32 @@ def maw_from_df( return maw -def _add_cellid(df, ds=None, gwf=None, x=None, y=None): +def _add_cellid(df, ds=None, gwf=None, x="x", y="y"): + """ + Intersect a DataFrame of point Data with the model grid, and add cellid-column. + + Parameters + ---------- + df : pd.DataFrame or gpd.GeoDataFrame + A (Geo)DataFrame containing the properties of the wells. + ds : xarray.Dataset + Dataset with model data. Either supply ds or gwf. The default is None. + gwf : flopy ModflowGwf + Groundwaterflow object. Only used when ds is None. The default is None. + x : str, optional + The column in df that contains the x-coordinate of the well. Only used when df + is a DataFrame. The default is 'x'. + y : str, optional + The column in df that contains the y-coordinate of the well. Only used when df + is a DataFrame. The default is 'y'. + + Returns + ------- + df : gpd.GeoDataFrame + A GeoDataFrame with a column named cellid that contains the icell2d-number + (vertex-grid) or (row, column) (structured grid). + + """ if not isinstance(df, gpd.GeoDataFrame): df = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df[x], df[y])) if "cellid" not in df.columns: @@ -192,17 +311,44 @@ def _add_cellid(df, ds=None, gwf=None, x=None, y=None): def _get_layer_multiplier_for_wells(df, top, botm, ds=None, gwf=None): + """ + Get factors (pandas.DataFrame) for each layer that all well screens intersects with. + + Parameters + ---------- + df : pd.DataFrame + A DataFrame containing the properties of the wells. + top : str + The column in df that contains the z-coordinate of the top of the well screen. + botm : str + The column in df that contains the z-coordinate of the bottom of the well + screen. + ds : xarray.Dataset + Dataset with model data. Either supply ds or gwf. The default is None. + gwf : flopy ModflowGwf + Groundwaterflow object. Only used when ds is None. The default is None. + + Returns + ------- + multipliers : pd.DataFrame + A DataFrame containg the multiplication factors, with the layers as the index + and the name of the well screens (the index of df) as columns. + + """ + """Get the factors """ # get required data either from gwf or ds - if ds is None: + if ds is not None: + ml_top = ds["top"].data + ml_bot = ds["botm"].data + kh = ds["kh"].data + layer = ds.layer + elif gwf is not None: ml_top = gwf.dis.top.array ml_bot = gwf.dis.botm.array kh = gwf.npf.k.array layer = range(gwf.dis.nlay.array) else: - ml_top = ds["top"].data - ml_bot = ds["botm"].data - kh = ds["kh"].data - layer = ds.layer + raise (Exception("Either supply ds or gwf to determine layer multiplyer")) multipliers = {} for index, irow in df.iterrows(): @@ -217,7 +363,30 @@ def _get_layer_multiplier_for_wells(df, top, botm, ds=None, gwf=None): def _get_layer_multiplier_for_well(cid, well_top, well_bot, ml_top, ml_bot, ml_kh): - """Get a factor for each layer that a well""" + """ + Get a factor (numpy array) for each layer that a well screen intersects with. + + Parameters + ---------- + cid : int or tuple of 2 ints + THe cellid of the well (either icell2d or (row, column). + well_top : float + The z-coordinate of the top of the well screen. + well_bot : float + The z-coordinate of the top of the well screen. + ml_top : numpy array + The top of the upper layer of the model (1d or 2d) + ml_bot : numpy array + The bottom of all cells of the model (2d or 3d) + ml_kh : numpy array + The horizontal conductivity of all cells of the model (2d or 3d). + + Returns + ------- + multiplier : numpy array + An array with a factor (between 0 and 1) for each of the model layers. + + """ # keep the tops and botms of the cell where the well is in ml_top_cid = ml_top[cid].copy() if isinstance(cid, int): @@ -244,6 +413,6 @@ def _get_layer_multiplier_for_well(cid, well_top, well_bot, ml_top, ml_bot, ml_k # the well does not cross any of the layers. Just return an array of zeros. multiplier = kd else: - # devide by the total kd to get a factor + # divide by the total kd to get a factor multiplier = kd / kd.sum() return multiplier From 72c6ae730af2f8930e3b7a04ef3ec59e8ce8ac43 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Thu, 5 Oct 2023 20:52:37 +0200 Subject: [PATCH 11/11] Solve codacy stuff and remove some leftover code --- nlmod/gwf/wells.py | 62 +++++++++++++++++++--------------------------- 1 file changed, 25 insertions(+), 37 deletions(-) diff --git a/nlmod/gwf/wells.py b/nlmod/gwf/wells.py index d9db0b81..87b556bf 100644 --- a/nlmod/gwf/wells.py +++ b/nlmod/gwf/wells.py @@ -46,8 +46,10 @@ def wel_from_df( The column in df that contains the z-coordinate of the bottom of the well screen. The defaults is 'botm'. Q : str, optional - The column in df containing the well discharge. This column can contain floats, - or strings belonging to timeseries added later. The default is "Q". + The column in df that contains the volumetric well rate. This column can contain + floats, or strings belonging to timeseries added later. A positive value + indicates recharge (injection) and a negative value indicates discharge + (extraction) The default is "Q". aux : str of list of str, optional The column(s) in df that contain auxiliary variables. The default is None. boundnames : str, optional @@ -160,8 +162,10 @@ def maw_from_df( The column in df that contains the z-coordinate of the bottom of the well screen. The defaults is 'botm'. Q : str, optional - The column in df containing the well discharge. This column can contain floats, - or strings belonging to timeseries added later. The default is "Q". + The column in df that contains the volumetric well rate. This column can contain + floats, or strings belonging to timeseries added later. A positive value + indicates recharge (injection) and a negative value indicates discharge + (extraction) The default is "Q". rw : str, optional The column in df that contains the radius for the multi-aquifer well. The default is "rw". @@ -194,40 +198,24 @@ def maw_from_df( df = _add_cellid(df, ds=ds, gwf=gwf, x=x, y=y) multipliers = _get_layer_multiplier_for_wells(df, top, botm, ds=ds, gwf=gwf) - maw_pakdata = [] - maw_conndata = [] - maw_perdata = [] + packagedata = [] + connectiondata = [] + perioddata = [] iw = 0 for index, irow in df.iterrows(): - try: - cid1 = gwf.modelgrid.intersect(irow[x], irow[y], irow[top], forgive=False) - cid2 = gwf.modelgrid.intersect(irow[x], irow[y], irow[botm], forgive=False) - except Exception: - logger.warning( - f"Well {index} outside of model domain ({irow[x]}, {irow[y]})" - ) - continue - kb = cid2[0] - if len(cid1) == 2: - kt, icell2d = cid1 - idomain_mask = gwf.modelgrid.idomain[kt : kb + 1, icell2d] > 0 - elif len(cid1) == 3: - kt, i, j = cid1 - idomain_mask = gwf.modelgrid.idomain[kt : kb + 1, i, j] > 0 - - wlayers = np.arange(kt, kb + 1)[idomain_mask] - wlayers = np.where(multipliers[index] > 0)[0] - # + + # [wellno, radius, bottom, strt, condeqn, ngwfnodes] pakdata = [iw, irow[rw], irow[botm], strt, condeqn, len(wlayers)] for iaux in aux: pakdata.append(irow[iaux]) if boundnames is not None: pakdata.append(irow[boundnames]) - maw_pakdata.append(pakdata) - # - maw_perdata.append([iw, "RATE", irow[Q]]) + packagedata.append(pakdata) + + # [wellno mawsetting] + perioddata.append([iw, "RATE", irow[Q]]) for iwellpart, k in enumerate(wlayers): if k == 0: @@ -249,7 +237,8 @@ def maw_from_df( scrn_top = np.min([irow[top], laytop]) scrn_bot = np.max([irow[botm], laybot]) - mawdata = [ + # [wellno, icon, cellid, scrn_top, scrn_bot, hk_skin, radius_skin] + condata = [ iw, iwellpart, cellid, @@ -258,7 +247,7 @@ def maw_from_df( 0.0, 0.0, ] - maw_conndata.append(mawdata) + connectiondata.append(condata) iw += 1 if len(aux) == 0: @@ -268,9 +257,9 @@ def maw_from_df( nmawwells=iw, auxiliary=aux, boundnames=boundnames is not None, - packagedata=maw_pakdata, - connectiondata=maw_conndata, - perioddata=maw_perdata, + packagedata=packagedata, + connectiondata=connectiondata, + perioddata=perioddata, **kwargs, ) @@ -312,7 +301,7 @@ def _add_cellid(df, ds=None, gwf=None, x="x", y="y"): def _get_layer_multiplier_for_wells(df, top, botm, ds=None, gwf=None): """ - Get factors (pandas.DataFrame) for each layer that all well screens intersects with. + Get factors (pandas.DataFrame) for each layer that well screens intersects with. Parameters ---------- @@ -335,7 +324,6 @@ def _get_layer_multiplier_for_wells(df, top, botm, ds=None, gwf=None): and the name of the well screens (the index of df) as columns. """ - """Get the factors """ # get required data either from gwf or ds if ds is not None: ml_top = ds["top"].data @@ -348,7 +336,7 @@ def _get_layer_multiplier_for_wells(df, top, botm, ds=None, gwf=None): kh = gwf.npf.k.array layer = range(gwf.dis.nlay.array) else: - raise (Exception("Either supply ds or gwf to determine layer multiplyer")) + raise (TypeError("Either supply ds or gwf to determine layer multiplyer")) multipliers = {} for index, irow in df.iterrows():