{ "cells": [ { "cell_type": "markdown", "id": "720d3f15", "metadata": {}, "source": [ "# The Lake Package\n", "The Lake (LAK) package can be used to model large surface water bodies with a single, uniform stage. The lake water balance is computed, and the resulting lake stage is determined by all inflows to and outflows from the lake. Several management options are available. Users can specify inflows and outflows, redirect discharge from other packages (e.g. DRN) to the lake, impose a fixed stage, or add a weir to control maximum water levels.\n", "\n", "In this example, the Henschotermeer lake in the Utrechtse Heuvelrug is modelled using the Lake package. For demonstration purposes, a weir is included with an invert level of 7.0 m MSL." ] }, { "cell_type": "code", "execution_count": null, "id": "a38194d9", "metadata": {}, "outputs": [], "source": [ "import os\n", "import pandas as pd\n", "import xarray as xr\n", "import matplotlib.pyplot as plt\n", "import matplotlib\n", "import flopy\n", "import nlmod" ] }, { "cell_type": "code", "execution_count": null, "id": "0d78cca8", "metadata": {}, "outputs": [], "source": [ "# set up pretty logging and show package versions\n", "nlmod.util.get_color_logger()\n", "nlmod.show_versions()" ] }, { "cell_type": "markdown", "id": "8ed27c72", "metadata": {}, "source": [ "## Define model extent and model workspace" ] }, { "cell_type": "code", "execution_count": null, "id": "86a246b9", "metadata": {}, "outputs": [], "source": [ "extent = [152_000, 156_000, 453_500, 456_000]\n", "model_name = \"henschotermeer\"\n", "model_ws = \"01_lake\"\n", "figdir, cachedir = nlmod.util.get_model_dirs(model_ws)" ] }, { "cell_type": "markdown", "id": "30cf20ed", "metadata": {}, "source": [ "## Download and prepare data\n", "We download:\n", "\n", "- the layers from the subsurface model Regis\n", "- the surface level data from AHN5\n", "- the layer 'waterdeel' from the Basisregistratie Grootschalige Toppgrafie (BGT)\n", "- level areas from waterboard Vallei en Veluwe" ] }, { "cell_type": "code", "execution_count": null, "id": "a4b40a04", "metadata": {}, "outputs": [], "source": [ "regis = nlmod.read.download_regis(extent, cachedir=cachedir, cachename=\"regis\")\n", "ahn = nlmod.read.ahn.download_ahn5(extent, cachedir=cachedir, cachename=\"ahn5\")\n", "bgt = nlmod.read.bgt.download_bgt(extent, cachedir=cachedir, cachename=\"bgt\")\n", "bgt = bgt.set_index(\"identificatie\")\n", "la = nlmod.read.waterboard.download_data(\n", " \"Vallei en Veluwe\",\n", " \"level_areas\",\n", " extent=extent,\n", " cachedir=cachedir,\n", " cachename=\"la_v_en_v\",\n", ")" ] }, { "cell_type": "markdown", "id": "29fc54d5", "metadata": {}, "source": [ "### Combine surface water data\n", "We add information from AHN and level areas to BGT-data." ] }, { "cell_type": "code", "execution_count": null, "id": "6e90397a", "metadata": {}, "outputs": [], "source": [ "columns = [\"summer_stage\", \"winter_stage\"]\n", "bgt = nlmod.util.gdf_intersection_join(\n", " la,\n", " bgt,\n", " columns=columns,\n", " min_total_overlap=0.0,\n", " desc=f\"Adding {columns} to bgt-data\",\n", ")\n", "\n", "bgt = nlmod.gwf.surface_water.add_min_ahn_to_gdf(bgt, ahn, buffer=5.0)\n", "bgt[\"summer_stage\"] = bgt[[\"summer_stage\", \"ahn_min\"]].max(1)\n", "bgt[\"winter_stage\"] = bgt[[\"winter_stage\", \"ahn_min\"]].max(1)\n", "bgt[\"stage\"] = bgt[[\"summer_stage\", \"winter_stage\"]].mean(1)\n", "\n", "bgt[\"lake\"] = None\n", "henschotermeer = bgt.area.idxmax() # 'W0662.ce4627405cb242ee8533946265644bb3'\n", "bgt.loc[henschotermeer, \"lake\"] = \"Henschotermeer\"" ] }, { "cell_type": "markdown", "id": "79ac6173", "metadata": {}, "source": [ "## Generate a model dataset\n", "We generate a model dataset with the same resolution as REGIS (100 x 100 m) and refine the model with level 2 (to 25x25 m) around the Henschotermeer." ] }, { "cell_type": "code", "execution_count": null, "id": "49981fdc", "metadata": {}, "outputs": [], "source": [ "ds = nlmod.to_model_ds(regis, model_name=model_name, model_ws=model_ws)\n", "\n", "# drop layers below \"PZWAz3\"\n", "ds = ds.sel(layer=ds.layer.loc[:\"PZWAz3\"])\n", "\n", "# refine around the edge of the lakes\n", "bgt_lake_boundary = bgt.loc[~bgt[\"lake\"].isna()].copy()\n", "bgt_lake_boundary.geometry = bgt_lake_boundary.boundary\n", "ds = nlmod.grid.refine(ds, refinement_features=[(bgt_lake_boundary, 2)])" ] }, { "cell_type": "markdown", "id": "9a144aa4", "metadata": {}, "source": [ "## Split surface water shapes by grid\n", "We split the surface water shapes by the modelgrid using the method `nlmod.grid.gdf_to_grid`. We then seperate the resulting features to a variable `lak_grid` for the LAK-package and a variable `drn_grid` for the DRN-pacakge." ] }, { "cell_type": "code", "execution_count": null, "id": "4c7e934c", "metadata": {}, "outputs": [], "source": [ "bgt_grid = nlmod.grid.gdf_to_grid(bgt, ds).set_index(\"cellid\")\n", "mask = bgt_grid[\"lake\"].isna()\n", "drn_grid = bgt_grid[mask].copy()\n", "lak_grid = bgt_grid[~mask].copy()\n", "assert not lak_grid.index.duplicated().any()\n", "# only keep cells where ahlf of the are of the cell is covered by lakes\n", "mask = lak_grid.area > 0.5 * ds[\"area\"].sel(icell2d=lak_grid.index)\n", "lak_grid = lak_grid[mask]\n", "# set the geometry to the entire cell\n", "gi = flopy.utils.GridIntersect(nlmod.grid.modelgrid_from_ds(ds))\n", "lak_grid.geometry = gi.geoms[lak_grid.index]\n", "\n", "# remove drains that overlap with the lake\n", "drn_grid = drn_grid.loc[~drn_grid.index.isin(lak_grid.index)]" ] }, { "cell_type": "markdown", "id": "f33c06df", "metadata": {}, "source": [ "## Improve model dataset\n", "### Set the top of the model\n", "The lake’s bottom elevation is defined by the top of the model. For this reason, we set the top of the model to 3.0 m MSL at the lake cells, which corresponds to the estimated lake bottom. The bottom elevation is particularly important in areas where the lake may dry out, which can also lead to convergence issues in the model. In this case, a bottom elevation of 3.0 m MSL is sufficiently deep, ensuring that the lake remains saturated throughout the transient simulation." ] }, { "cell_type": "code", "execution_count": null, "id": "ea5adb13", "metadata": {}, "outputs": [], "source": [ "top = nlmod.resample.structured_da_to_ds(nlmod.resample.fillnan_da(ahn), ds)\n", "# The botom of the lake is at 3.0 m NAP\n", "top[lak_grid.index] = 3.0\n", "ds = nlmod.layers.set_model_top(ds, top)" ] }, { "cell_type": "markdown", "id": "2b8053e0", "metadata": {}, "source": [ "### Set time dimension\n", "Before downloading meteorological data, we define the time dimension of our model dataset. The simulation covers five years, from the beginning of 2020 to the beginning of 2025, using monthly stress periods. The model is initialized with a steady-state stress period representing the mean meteorological conditions of 2019." ] }, { "cell_type": "code", "execution_count": null, "id": "ea6d113a", "metadata": {}, "outputs": [], "source": [ "time = pd.date_range(\"2020\", \"2025\", freq=\"MS\")\n", "ds = nlmod.time.set_ds_time(ds, start=\"2019\", time=time)" ] }, { "cell_type": "markdown", "id": "3b5d66d4", "metadata": {}, "source": [ "### Download recharge data\n", "We use the method `nlmod.read.knmi.get_recharge` to download meteorological data. Instead of calculating recharge as precipitation minus evaporation, we retrieve both variables separately by setting `method=\"separate\"`. Since our model area is relatively small, we represent it using a single meteorological and precipitation station (`most_common_station=True`). The variables recharge and evaporation are defined per station rather than for each model cell, which is achieved with `add_stn_dimensions=True`. Finally, to adopt the new default and suppress warnings about hourly precision, we set `hourly_precision=True`." ] }, { "cell_type": "code", "execution_count": null, "id": "29762f4a", "metadata": {}, "outputs": [], "source": [ "rch_ds = nlmod.read.knmi.get_recharge(\n", " ds=ds,\n", " method=\"separate\",\n", " most_common_station=True,\n", " add_stn_dimensions=True,\n", " hourly_precision=True,\n", " cachedir=cachedir,\n", " cachename=\"knmi\",\n", ")\n", "ds.update(rch_ds)\n", "\n", "ds[\"starting_head\"] = xr.full_like(ds.botm, 5.0)" ] }, { "cell_type": "markdown", "id": "018c15f1", "metadata": {}, "source": [ "## Generate and run model\n", "### Generate a FloPy sim and gwf " ] }, { "cell_type": "code", "execution_count": null, "id": "4e28c259", "metadata": {}, "outputs": [], "source": [ "sim = nlmod.sim.sim(ds) # simulation\n", "tdis = nlmod.sim.tdis(ds, sim) # time discretization\n", "ims = nlmod.sim.ims(\n", " sim,\n", " complexity=\"COMPLEX\",\n", " inner_dvclose=0.01, # needed so lake balance is correct\n", " outer_dvclose=0.01, # needed so lake balance is correct\n", " #rcloserecord=[0.01, \"STRICT\"],\n", ") # ims solver\n", "gwf = nlmod.gwf.gwf(ds, sim, under_relaxation=True) # groundwater flow model\n", "dis = nlmod.gwf.dis(ds, gwf) # spatial discretization\n", "npf = nlmod.gwf.npf(ds, gwf) # node property flow\n", "sto = nlmod.gwf.sto(ds, gwf) # storage\n", "ic = nlmod.gwf.ic(ds, gwf) # initial conditions\n", "oc = nlmod.gwf.oc(ds, gwf) # output control" ] }, { "cell_type": "markdown", "id": "01f96640", "metadata": {}, "source": [ "### Meteorological stresses" ] }, { "cell_type": "code", "execution_count": null, "id": "017a1f89", "metadata": {}, "outputs": [], "source": [ "rch = nlmod.gwf.rch(ds, gwf)\n", "evt = nlmod.gwf.evt(ds, gwf)" ] }, { "cell_type": "markdown", "id": "95a60507", "metadata": {}, "source": [ "### Lake package for the Henschotermeer" ] }, { "cell_type": "code", "execution_count": null, "id": "d9ba825c", "metadata": {}, "outputs": [], "source": [ "# %% \n", "lak_resistance = 10.0 # days\n", "lak_grid[\"clake\"] = lak_resistance\n", "lak_grid[\"strt\"] = 6.0\n", "lak_grid[\"lakeout\"] = -1\n", "lak_grid[\"outlet_invert\"] = 7.0\n", "nlmod.gwf.lake_from_gdf(\n", " gwf,\n", " lak_grid,\n", " ds,\n", " rainfall=ds[\"recharge\"].to_pandas().iloc[:, 0],\n", " evaporation=ds[\"evaporation\"].to_pandas().iloc[:, 0] * 1.25,\n", " boundname_column=\"lake\",\n", " #maximum_stage_change=0.000001,\n", " #maximum_iterations=200,\n", ");" ] }, { "cell_type": "markdown", "id": "8eab2e92", "metadata": {}, "source": [ "### Drain package for all other surface water" ] }, { "cell_type": "code", "execution_count": null, "id": "85201dec", "metadata": {}, "outputs": [], "source": [ "drn_resistance = 1.0 # days\n", "drn_grid[\"cond\"] = drn_grid.area / drn_resistance\n", "spd = nlmod.gwf.surface_water.build_spd(drn_grid, \"DRN\", ds)\n", "drn = flopy.mf6.ModflowGwfdrn(gwf, stress_period_data={0: spd})" ] }, { "cell_type": "markdown", "id": "6eeb21c1", "metadata": {}, "source": [ "### Write input-files and run Modflow 6" ] }, { "cell_type": "code", "execution_count": null, "id": "c206be01", "metadata": {}, "outputs": [], "source": [ "nlmod.sim.write_and_run(sim, ds)" ] }, { "cell_type": "markdown", "id": "42a7f0e1", "metadata": {}, "source": [ "## Post-processing\n", "### Plot the head in the top layer\n", "We plot the time-averaged head in the first model layer, `BXz2`. This layer is absent in the western part of the model domain, which is why this area appears white in the figure. Surface water features (`bgt`) are plotted on top of the map for reference. The finer model cells along the shoreline of Lake Henschotermeer are clearly visible." ] }, { "cell_type": "code", "execution_count": null, "id": "abbfb66b", "metadata": {}, "outputs": [], "source": [ "head = nlmod.gwf.get_heads_da(ds)\n", "ax = nlmod.plot.map_array(head.sel(layer=\"BXz2\").mean(\"time\"), ds=ds)\n", "bgt.plot(edgecolor='k', facecolor='none', ax=ax);" ] }, { "cell_type": "markdown", "id": "cfead887", "metadata": {}, "source": [ "### Plot the Lake stage\n", "If the parameter `boundname_column` is specified in `nlmod.gwf.lake_from_gdf`, a CSV file named `lak_STAGE.csv` is generated containing observations of lake stages. Below, we plot the lake stage. Due to the wet conditions during the winter of 2023/2024, the lake level rises to 7.0 m MSL. At this elevation, the weir becomes active, preventing any further significant increase in lake stage." ] }, { "cell_type": "code", "execution_count": null, "id": "64b8e5db", "metadata": {}, "outputs": [], "source": [ "lak_stage = pd.read_csv(os.path.join(ds.model_ws, \"lak_STAGE.csv\"), index_col=0)\n", "lak_stage.index = pd.to_datetime(ds.time.start) + pd.to_timedelta(lak_stage.index, \"d\")\n", "lak_stage.columns = [x.capitalize() for x in lak_stage.columns]\n", "\n", "f, ax = plt.subplots(figsize=(10, 6), layout=\"constrained\")\n", "lak_stage.plot(ax=ax)\n", "ax.set_xlabel(\"\")\n", "ax.set_ylabel(\"lake stage (m MSL)\")\n", "ax.grid()" ] }, { "cell_type": "markdown", "id": "9f3f7169", "metadata": {}, "source": [ "### Plot a cross-section\n", "We plot a west–east cross-section at y = 454 612, showing the model layers. We can see the Utrechtse Heuvelrug on the left (west), with decreasing surface level eastwards. The lake \"Henschotermeer\" is displayed in blue, using the mean lake stage over the simulation period." ] }, { "cell_type": "code", "execution_count": null, "id": "77a9d941", "metadata": {}, "outputs": [], "source": [ "f, ax = plt.subplots(figsize=(10, 6), layout=\"constrained\")\n", "line = [(extent[0], 454_612), (extent[1], 454_612)]\n", "dcs = nlmod.plot.DatasetCrossSection(ds, line, ax=ax, zmin=-100)\n", "dcs.plot_layers(colors=nlmod.read.regis.get_legend())\n", "dcs.plot_grid(vertical=False, lw=0.5)\n", "dcs.label_layers()\n", "\n", "lak_gdf = lak_grid.dissolve(\"lake\")\n", "\n", "\n", "def plot_lake_in_cs(linestring, dcs, height, ax, color=\"C0\", zorder=0):\n", " xy = (dcs.line.project(linestring.boundary.geoms[0]), dcs.zmin)\n", " width = dcs.line.project(linestring.boundary.geoms[1]) - xy[0]\n", " rect = matplotlib.patches.Rectangle(\n", " xy, width, height, color=color, zorder=zorder\n", " )\n", " ax.add_patch(rect)\n", "\n", "# plot the lake\n", "for lake in bgt.loc[~bgt.lake.isna(), \"lake\"].unique():\n", " mean_lake_stage = lak_stage[lake].mean()\n", " ylim = ax.get_ylim()\n", " height = mean_lake_stage - ylim[0]\n", " mask = bgt.lake == lake\n", " for geom in lak_gdf.loc[[lake]].intersection(dcs.line).geometry.values:\n", " if geom.geom_type == \"MultiLineString\":\n", " for geom_line in geom.geoms:\n", " plot_lake_in_cs(geom_line, dcs, height, ax)\n", " else:\n", " plot_lake_in_cs(geom, dcs, height, ax)\n", "\n", "\n", "axes_bounds = nlmod.plot.get_inset_map_bounds(ax, extent, height=0.3)\n", "mapax = nlmod.plot.inset_map(ax, extent, axes_bounds=axes_bounds)\n", "nlmod.plot.add_xsec_line_and_labels(line, ax, mapax)\n", "ax.set_xlabel(\"x (length along cross-section, m)\")\n", "ax.set_ylabel(\"z (m MSL)\");" ] }, { "cell_type": "markdown", "id": "4eae01d9", "metadata": {}, "source": [ "### Plot a water balance of the lake\n", "The main advantage of using the LAK package is that the lake water balance is explicitly accounted for. This water balance is saved in the file `lak.bgt`. Below, we read this file and plot the results. The primary fluxes into and out of the lake are rainfall (blue) and evaporation (orange), which are model inputs. The calculated flux between the lake and the groundwater (GWF, brown) shows a substantial contribution of groundwater flow to the lake during the winter of 2023/2024. This inflow, and rainfall, causes the lake stage to rise, resulting in a negative storage term (grey). From the winter of 2023/2024 onward, the weir becomes active, causing external outflow of water (red)." ] }, { "cell_type": "code", "execution_count": null, "id": "807e26aa", "metadata": {}, "outputs": [], "source": [ "fname = os.path.join(ds.model_ws, \"lak.bgt\")\n", "cbf_lak = nlmod.gwf.output.get_cellbudgetfile(fname=fname, ds=ds)\n", "\n", "# read budgets\n", "lake_bgt = {}\n", "for rn in cbf_lak.get_unique_record_names():\n", " key = rn.strip().decode()\n", " lake_bgt[key] = cbf_lak.get_data(text=rn)\n", "\n", "colors = {\n", " \"GWF\": \"tab:brown\",\n", " \"RAINFALL\": \"tab:blue\",\n", " \"CONSTANT\": \"tab:purple\",\n", " \"EVAPORATION\": \"tab:orange\",\n", " \"STORAGE\": \"tab:gray\",\n", " \"FROM-MVR\": \"tab:green\",\n", " \"EXT-OUTFLOW\": \"tab:red\",\n", "}\n", "\n", "\n", "for name in lak_grid[\"lake\"].unique():\n", " lakeno = int(lak_grid.loc[lak_grid[\"lake\"] == name, \"lakeno\"].iloc[0])\n", " # generate a DataFrame\n", " data = {}\n", " for key in lake_bgt.keys():\n", " data[key] = [x[x[\"node\"] == lakeno + 1][\"q\"].sum() for x in lake_bgt[key]]\n", "\n", " index = pd.to_datetime(ds.time.start) + pd.to_timedelta(cbf_lak.get_times(), \"d\")\n", " df = pd.DataFrame(data, index=index)\n", "\n", " f, ax = plt.subplots(figsize=(10, 6), layout=\"constrained\")\n", "\n", " dfs = df.copy()\n", " dfs.index = [pd.to_datetime(ds.time.start)] + list(df.index[:-1])\n", " df_block = pd.concat((df, dfs)).sort_index(kind=\"mergesort\")\n", "\n", " inflow = df_block.where(df_block > 0, 0.0)\n", " inflow = inflow.loc[:, ~(inflow == 0).all()]\n", "\n", " outflow = df_block.where(df_block < 0, 0.0)\n", " outflow = outflow.loc[:, ~(outflow == 0).all()]\n", "\n", " ymax = max(inflow.sum(1).max(), outflow.sum(1).max())\n", "\n", " inflow.plot.area(color=colors, ax=ax, linewidth=0)\n", " handles_in, labels_in = ax.get_legend_handles_labels()\n", " ax.set_ylim(-ymax, ymax)\n", " outflow.plot.area(color=colors, ax=ax, linewidth=0)\n", " ax.set_ylim(-ymax, ymax)\n", "\n", " ax.set_xlim(ds.time.data[0], ds.time.data[-1])\n", " ax.set_ylabel(\"Outflow (negative) and inflow (positive), in m3/d\")\n", " nlmod.plot.title_inside(name, ax=ax)\n", "\n", " # remove double legend entries\n", " handles_all, labels_all = ax.get_legend_handles_labels()\n", "\n", " handles = handles_in[::-1]\n", " labels = labels_in[::-1]\n", "\n", " # move storage to the middle of the legend\n", " index = labels.index(\"STORAGE\")\n", " handles.append(handles.pop(index))\n", " labels.append(labels.pop(index))\n", " for handle, label in zip(handles_all, labels_all):\n", " if label not in labels:\n", " handles.append(handle)\n", " labels.append(label)\n", "\n", " ax.legend(handles, labels, loc=2)\n", " ax.grid()" ] } ], "metadata": { "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 5 }