<jupyter_start><jupyter_text>NLmod examples  This example notebook shows a basic example of a model created with the `nlmod` package. The `nlmod` package contains functions to create modflow models anywhere in the Netherlands.  Contents1. [Create model](create)2. [Run model](run)3. [Visualise](visualise)<jupyter_code>import matplotlib.pyplot as plt
import flopy
import os
import geopandas as gpd

import nlmod

import logging
print(f'nlmod version: {nlmod.__version__}')

# toon informatie bij het aanroepen van functies
logging.basicConfig(level=logging.INFO)<jupyter_output>nlmod version: 0.2.0<jupyter_text>[1. Create model](TOC)With the code below we create a modflow model with the name 'IJmuiden'. This model has the following properties :- an extent that covers part of the Northsea, Noordzeekanaal and the small port city IJmuiden.- a structured grid based on the subsurface models [Regis](https://www.dinoloket.nl/regis-ii-het-hydrogeologische-model) and [Geotop](https://www.dinoloket.nl/detaillering-van-de-bovenste-lagen-met-geotop). The Regis layers that are not present within the extent are removed. In this case we use b'MSz1' as the bottom layer of the model. Use `nlmod.regis.get_layer_names()` to get all the layer names of Regis. All Regis layers below this layer are not used in the model. Geotop is used to replace the holoceen layer in Regis because there is no kh or kv defined for the holoceen in Regis. Part of the model is in the North sea. Regis and Geotop have no data there. Therefore the Regis and Geotop layers are extrapolated from the shore and the seabed is added using bathymetry data from [Jarkus](https://www.openearth.nl/rws-bathymetry/2018.html).- starting heads of 1 in every cell.- the model is a steady state model of a single time step.- big surface water bodies (Northsea, IJsselmeer, Markermeer, Noordzeekanaal) within the extent are added as a general head boundary. The surface water bodies are obtained from a [shapefile](..\data\opp_water.shp).- surface drainage is added using [ahn](https://www.ahn.nl) data and a default conductance of $1000 m^2/d$- recharge is added using data from the [knmi](https://www.knmi.nl/nederland-nu/klimatologie/daggegevens) using the following steps:~~    1. Check for each cell which KNMI weather and/or rainfall station is closest.    2. Download the data for the stations found in 1. for the model period. For a steady state stress period the average precipitation and evaporation of 8 years before the stress period time is used.    3. Combine precipitation and evaporation data from step 2 to create a recharge time series for each cell    4. Add the timeseries to the model dataset and create the recharge package.- constant head boundaries are added to the model edges in every layer. The starting head is used as constant head.<jupyter_code># model settings
model_ws = 'model1'
model_name = 'IJmuiden'
use_cache = True
extent = [95000., 105000., 494000., 500000.]
delr = 100.
delc = 100.
steady_state = True
start_time = '2015-1-1'
gridtype = 'structured'
use_regis = True
regis_botm_layer = b'MSz1'
use_geotop = True
add_northsea = True
starting_head = 1.0
# create empty model dataset
model_ds = nlmod.mdims.get_empty_model_ds(model_name, model_ws)

# add time discretisation
model_ds = nlmod.mdims.set_model_ds_time(model_ds,
                                         start_time,
                                         steady_state, perlen=365*5)

# add spatial discretisation
extent, nrow, ncol = nlmod.read.regis.fit_extent_to_regis(extent,
                                                          delr,
                                                          delc)

# layer model
layer_model = nlmod.read.regis.get_combined_layer_models(extent, delr, delc,
                                                         use_regis=use_regis,
                                                         regis_botm_layer=regis_botm_layer,
                                                         use_geotop=use_geotop,
                                                         cachedir=model_ds.cachedir,
                                                         cachename='combined_layer_ds.nc')

# create modflow packages
sim, gwf = nlmod.mfpackages.sim_tdis_gwf_ims_from_model_ds(model_ds)
# update model_ds from layer model
model_ds = nlmod.mdims.update_model_ds_from_ml_layer_ds(model_ds,
                                                        layer_model,
                                                        keep_vars=['x', 'y'],
                                                        gridtype=gridtype,
                                                        add_northsea=add_northsea,
                                                        cachedir=model_ds.cachedir)

# Create discretization
nlmod.mfpackages.dis_from_model_ds(model_ds, gwf)

# create node property flow
nlmod.mfpackages.npf_from_model_ds(model_ds, gwf)

# Create the initial conditions package
nlmod.mfpackages.ic_from_model_ds(model_ds, gwf, starting_head=starting_head)

# Create the output control package
nlmod.mfpackages.oc_from_model_ds(model_ds, gwf);
# voeg grote oppervlaktewaterlichamen toe o.b.v. rws shape
da_name = 'rws_oppwater'
rws_ds = nlmod.read.rws.get_surface_water(model_ds,
                                           da_name,
                                           cachedir=model_ds.cachedir,
                                           cachename=da_name)
model_ds.update(rws_ds)
ghb = nlmod.mfpackages.ghb_from_model_ds(model_ds, gwf, da_name)
# surface level drain
ahn_ds = nlmod.read.ahn.get_ahn(model_ds, cachedir=model_ds.cachedir, cachename='ahn')
model_ds.update(ahn_ds)

drn = nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf)


# add constant head cells at model boundaries
model_ds.update(nlmod.mfpackages.constant_head.get_chd_at_model_edge(model_ds, model_ds['idomain']))    
chd = nlmod.mfpackages.chd_from_model_ds(model_ds, gwf, head='starting_head')
# add knmi recharge to the model datasets
knmi_ds = nlmod.read.knmi.get_recharge(model_ds, cachedir=model_ds.cachedir, cachename='recharge')
model_ds.update(knmi_ds)

# create recharge package
rch = nlmod.mfpackages.rch_from_model_ds(model_ds, gwf)<jupyter_output>INFO:nlmod.cache:using cached data -> recharge.nc<jupyter_text>A big part of the model data is stored in the variable `model_ds` which is an `xarray.Dataset`. The data is shown below.<jupyter_code>model_ds<jupyter_output><empty_output><jupyter_text>[2. Write and Run](TOC)Now that we've created all the modflow packages we need to write them to modflow files. You always have to write the modflow data to the model workspace before you can run the model. You can write the model files and run the model using the function `nlmod.util.write_and_run_model()` as shown below. This function has two additional options:1. Write the model dataset to the disk if `write_model_ds` is `True`. This makes it easier and faster to load model data if you ever need it. 2. Write a copy of this Jupyter Notebook to the same directory as the modflow files if `nb_path` is the name of this Jupyter Notebook. It can be useful to have a copy of the script that created the modflow files, together with the files.<jupyter_code>nlmod.util.write_and_run_model(gwf, model_ds, write_model_ds=True, nb_path='01_basic_model.ipynb')<jupyter_output>INFO:nlmod.util:write script 2021121001_basic_model.ipynb to model workspace
INFO:nlmod.util:write model dataset to cache
INFO:nlmod.util:write modflow files to model workspace<jupyter_text>[3. Visualise](TOC)Using the `model_ds` and `gwf` variables it is quite easy to visualise model data. Below the modelgrid together with the surface water is shown.<jupyter_code>ax = nlmod.visualise.plots.plot_modelgrid(model_ds, gwf)
ax.figure.savefig(os.path.join(model_ds.figdir, 'mgrid_swater.png'), bbox_inches='tight')<jupyter_output>ERROR:fiona._env:PROJ: proj_identify: Cannot find proj.db
INFO:fiona.ogrext:Failed to auto identify EPSG: 7<jupyter_text>Data from a model with a structured grid can be easily visualised using the model dataset. Below some examples<jupyter_code>fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(14, 11))
model_ds['ahn'].plot(ax=axes[0][0])
model_ds['bot'][0].plot(ax=axes[0][1])
model_ds['idomain'][0].plot(ax=axes[1][0])
model_ds['chd'][0].plot(ax=axes[1][1])
for axes1 in axes:
    for ax in axes1:
        ax.axis('scaled')

fig.savefig(os.path.join(model_ds.figdir, 'ahn_bot_idom_chd.png'), bbox_inches='tight')

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(14, 11))
model_ds['bathymetry'].plot(ax=axes[0][0])
model_ds['northsea'].plot(ax=axes[0][1])
model_ds['kh'][1].plot(ax=axes[1][0])
model_ds['recharge'].plot(ax=axes[1][1])

for axes1 in axes:
    for ax in axes1:
        ax.axis('scaled')
fig.savefig(os.path.join(model_ds.figdir, 'bath_nsea_kh_top.png'), bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0844833389431168, 'fraction_numerical': 0.013575675978907215, 'mean_word_length': 3.333495381623724, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 22, 'https://': 5, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 35, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '663924', 'license': 'MIT', 'n_tokens_mistral': 2861, 'n_tokens_neox': 2744, 'n_words': 835, 'path': 'examples/01_basic_model.ipynb', 'repo_name': 'ArtesiaWater/nlmod', 'size': '8913'}