From 9135395e7c2ce9254075d4f1af186d2f2878802d Mon Sep 17 00:00:00 2001 From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com> Date: Wed, 30 Oct 2024 13:00:32 +0100 Subject: [PATCH] rewrap text --- docs/changelog.qmd | 12 +- docs/developments/guide.qmd | 7 +- docs/developments/julia_structs.qmd | 38 +- docs/getting_started/building_a_model.qmd | 96 ++- .../download_example_models.qmd | 10 +- docs/getting_started/install.qmd | 74 ++- docs/getting_started/running_wflow.qmd | 13 +- docs/home/case_studies.qmd | 101 ++-- docs/home/publications.qmd | 236 ++++---- docs/index.qmd | 36 +- docs/model_docs/index.qmd | 20 +- docs/model_docs/lateral/gwf.qmd | 182 +++--- docs/model_docs/lateral/kinwave.qmd | 129 ++-- docs/model_docs/lateral/local-inertial.qmd | 164 +++--- docs/model_docs/lateral/sediment_flux.qmd | 361 ++++++------ docs/model_docs/lateral/waterbodies.qmd | 74 +-- docs/model_docs/model_configurations.qmd | 92 +-- docs/model_docs/parameters_lateral.qmd | 201 ++++--- docs/model_docs/parameters_vertical.qmd | 91 ++- docs/model_docs/vertical/sbm.qmd | 549 +++++++++--------- docs/model_docs/vertical/sediment.qmd | 197 ++++--- docs/model_docs/vertical/shared_processes.qmd | 101 ++-- docs/user_guide/bmi.qmd | 37 +- docs/user_guide/faq.qmd | 11 +- docs/user_guide/fews.qmd | 19 +- docs/user_guide/multi_threading.qmd | 39 +- docs/user_guide/output.qmd | 4 +- docs/user_guide/required_files.qmd | 52 +- docs/user_guide/toml_file.qmd | 168 +++--- docs/user_guide/warm_states.qmd | 23 +- docs/user_guide/zmq_server.qmd | 4 +- 31 files changed, 1577 insertions(+), 1564 deletions(-) diff --git a/docs/changelog.qmd b/docs/changelog.qmd index 65a90d606..d9244485d 100644 --- a/docs/changelog.qmd +++ b/docs/changelog.qmd @@ -4,8 +4,8 @@ title: "Changelog" All notable changes to this project will be documented in this file. -The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), -and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this +project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## v1.0.0 @@ -21,7 +21,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## v0.x releases -### [unreleased] +### v0.8.1 - 2024-08-27 + +#### Fixed +- Reduce allocations in update of vertical `SBM` concept. + +### v0.8.0 - 2024-08-19 #### Fixed - Added missing BMI function `get_grid_size`, it is used for unstructured grids, for example @@ -115,6 +120,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Local inertial routing to `sbm_gwf` model type. - Water demand and allocation computations for model types `sbm` and `sbm_gwf`. + ### v0.7.3 - 2024-01-12 #### Fixed diff --git a/docs/developments/guide.qmd b/docs/developments/guide.qmd index 92bd93500..5f027754e 100644 --- a/docs/developments/guide.qmd +++ b/docs/developments/guide.qmd @@ -4,9 +4,10 @@ title: Developers guide ## Contributions and reporting issues -We welcome reporting of issues on [our GitHub page](https://github.com/Deltares/Wflow.jl/issues/new/choose). -Please provide a minimum working example so we are able to reproduce the issue. Furthermore, we -welcome contributions. We follow the [ColPrac guide for collaborative practices](https://github.com/SciML/ColPrac). New +We welcome reporting of issues on [our GitHub +page](https://github.com/Deltares/Wflow.jl/issues/new/choose). Please provide a minimum working +example so we are able to reproduce the issue. Furthermore, we welcome contributions. We follow +the [ColPrac guide for collaborative practices](https://github.com/SciML/ColPrac). New contributors should make sure to read that guide. ## Style/decisions diff --git a/docs/developments/julia_structs.qmd b/docs/developments/julia_structs.qmd index 3627e8b13..0ed155b9d 100644 --- a/docs/developments/julia_structs.qmd +++ b/docs/developments/julia_structs.qmd @@ -4,8 +4,8 @@ title: Julia structures ## Model - Below the composite type that represents all different aspects of a `Wflow.Model`, such as - the network, parameters, clock, model type, configuration and input and output. +Below the composite type that represents all different aspects of a `Wflow.Model`, such as the +network, parameters, clock, model type, configuration and input and output. ```julia struct Model{N,L,V,R,W,T} @@ -21,9 +21,9 @@ end ``` The `lateral` field of the `struct Model` can contain different lateral concepts. For each -wflow model these different lateral concepts are mapped through the use of a `NamedTuple`. -The `vertical` field of the `struct Model` always contains one vertical concept, for example -the SBM vertical concept. +wflow model these different lateral concepts are mapped through the use of a `NamedTuple`. The +`vertical` field of the `struct Model` always contains one vertical concept, for example the +SBM vertical concept. Below an example how lateral concepts are mapped for the SBM model through a `NamedTuple`: @@ -32,22 +32,21 @@ Below an example how lateral concepts are mapped for the SBM model through a `Na ``` The `subsurface` part is mapped to the lateral subsurface flow kinematic wave concept, the -`land` part is mapped the overland flow kinematic wave concept and the `river` part is -mapped to the river flow kinematic wave concept. Knowledge of this specific mapping is -required to understand and correctly set input, output and state variables in the TOML -configuration file, see also [Config and TOML](@ref config_toml). This mapping is described in more -detail for each model in the section Models. Also the `struct` of each mapped concept is -provided, so one can check the internal variables in the code. These structs are defined as -a parametric composite type, with type parameters between curly braces after the `struct` -name. See also the next paragraph [Vertical and lateral models](@ref) for a more -detailed description. +`land` part is mapped the overland flow kinematic wave concept and the `river` part is mapped +to the river flow kinematic wave concept. Knowledge of this specific mapping is required to +understand and correctly set input, output and state variables in the TOML configuration file, +see also [Config and TOML](@ref config_toml). This mapping is described in more detail for each +model in the section Models. Also the `struct` of each mapped concept is provided, so one can +check the internal variables in the code. These structs are defined as a parametric composite +type, with type parameters between curly braces after the `struct` name. See also the next +paragraph [Vertical and lateral models](@ref) for a more detailed description. ## Vertical and lateral models The different model concepts used in wflow are defined as parametric [composite types](https://docs.julialang.org/en/v1/manual/types/#Composite-Types). For example the vertical `SBM` concept is defined as follows: `struct SBM{T,N,M}`. `T`, `N` and `M` between -curly braces after the `struct` name refer to type parameters, for more information about -type parameters you can check out [Type +curly braces after the `struct` name refer to type parameters, for more information about type +parameters you can check out [Type parameters](https://docs.julialang.org/en/v1/manual/types/#man-parametric-composite-types). Since these parameters can be of any type, it is possible to declare an unlimited number of composite types. The type parameters are used to set the type of `struct` fields, below an @@ -78,10 +77,9 @@ example with a part of the `SBM` struct: ``` The type parameter `T` is used in wflow as a subtype of `AbstractFloat`, allowing to store -fields with a certain floating point precision (e.g. `Float64` or `Float32`) in a flexible -way. `N` refers to the maximum number of soil layers of the `SBM` soil column, and `M` -refers to the maximum number of soil layers + 1. See also part of the following instance of -`SBM`: +fields with a certain floating point precision (e.g. `Float64` or `Float32`) in a flexible way. +`N` refers to the maximum number of soil layers of the `SBM` soil column, and `M` refers to the +maximum number of soil layers + 1. See also part of the following instance of `SBM`: ```julia sbm = SBM{Float,maxlayers,maxlayers + 1}( diff --git a/docs/getting_started/building_a_model.qmd b/docs/getting_started/building_a_model.qmd index c094091b2..5f7b1c8b0 100644 --- a/docs/getting_started/building_a_model.qmd +++ b/docs/getting_started/building_a_model.qmd @@ -4,13 +4,12 @@ title: Building a model from scratch ## HydroMT-wflow -[hydroMT](https://github.com/Deltares/hydromt) is a Python package, developed by Deltares, -to build and analyze hydro models. It provides a generic model api with attributes to -access the model schematization, (dynamic) forcing data, results and states. +[hydroMT](https://github.com/Deltares/hydromt) is a Python package, developed by Deltares, to +build and analyze hydro models. It provides a generic model api with attributes to access the +model schematization, (dynamic) forcing data, results and states. -The wflow plugin -[hydroMT-wflow](https://github.com/Deltares/hydromt_wflow) of hydroMT can be used to build -and analyze the following model configurations: +The wflow plugin [hydroMT-wflow](https://github.com/Deltares/hydromt_wflow) of hydroMT can be +used to build and analyze the following model configurations: - [wflow\_sbm + kinematic wave routing](../model_docs/model_configurations.html#sbm-kinematic-wave) - [wflow\_sbm + local inertial river and floodplain](../model_docs/model_configurations.html#sbm-local-inertial-river) @@ -20,9 +19,9 @@ and analyze the following model configurations: To learn more about the wflow plugin of this Python package, we refer to the [hydroMT-wflow documentation](https://deltares.github.io/hydromt_wflow/latest/index.html). -To inspect or modify (for example in QGIS) the netCDF static data of these wflow models it -is convenient to export the maps to a raster format. This can be done as part of the -hydroMT-wflow plugin, see also the following [example] +To inspect or modify (for example in QGIS) the netCDF static data of these wflow models it is +convenient to export the maps to a raster format. This can be done as part of the hydroMT-wflow +plugin, see also the following [example] (https://deltares.github.io/hydromt_wflow/latest/_examples/convert_staticmaps_to_mapstack.html). It is also possible to create again the netCDF static data file based on the modified raster map stack. @@ -30,9 +29,9 @@ map stack. ## Data requirements The actual data requirements depend on the application of the Model and the Model type. Both -forcing and static data should be provided in netCDF format, with the same grid definition -for forcing and static data. The only exception is storage and rating curves for lakes, that -should be provided in CSV format, see also [Additional settings](@ref). +forcing and static data should be provided in netCDF format, with the same grid definition for +forcing and static data. The only exception is storage and rating curves for lakes, that should +be provided in CSV format, see also [Additional settings](@ref). * Forcing data: - Precipitation @@ -52,15 +51,14 @@ An approach to generate `ldd` data is to make use of the Python package [pyflwdir](https://github.com/Deltares/pyflwdir): + to [upscale existing flow direction - data](https://deltares.github.io/pyflwdir/latest/_examples/upscaling.html) as the 3 arcsec MERIT - Hydro data (Yamazaki et al., 2019) + data](https://deltares.github.io/pyflwdir/latest/_examples/upscaling.html) as the 3 arcsec + MERIT Hydro data (Yamazaki et al., 2019) + or to [derive flow directions from elevation data](https://deltares.github.io/pyflwdir/latest/_examples/from_dem.html), -see also Eilander et al. (2021) for more information. -Pyflwdir is also used by the [hydroMT](@ref) Python package described in the next paragraph. -Another approach to generate `ldd` data is to make use of PCRaster functionality, see for -example +see also Eilander et al. (2021) for more information. Pyflwdir is also used by the +[hydroMT](@ref) Python package described in the next paragraph. Another approach to generate +`ldd` data is to make use of PCRaster functionality, see for example [lddcreate](https://pcraster.geo.uu.nl/pcraster/4.3.1/documentation/pcraster_manual/sphinx/op_lddcreate.html). Optionally, but also not directly part of a model component are `gauge` locations, that are @@ -73,24 +71,24 @@ section of Model parameters. For wflow\_sbm models there are two ways to include flow: 1. The kinematic wave approach (see section [Subsurface flow routing](@ref)) as part of the - `sbm` model type. Parameters that are part of this component are described in the - [Lateral subsurface flow](@ref params_ssf) section of Model parameters. Input parameters - for this component are derived from the SBM vertical concept and the land slope. One - external parameter [`ksathorfrac`](@ref params_ssf) is used to calculate the horizontal - hydraulic conductivity at the soil surface `kh_0`. -2. Groundwater flow (see section [Groundwater flow component](@ref lateral_gwf)) as part of - the `sbm_gwf` model type. For the unconfined aquifer the input parameters are described - in the section [Unconfined aquifer](@ref) of Model parameters. The bottom (`bottom`) of - the groundwater layer is derived from from the `soilthickness` [mm] parameter of `SBM` - and the provided surface elevation `altitude` [m] as part of the static input. The `area` - parameter is derived from the model grid. Parameters that are part of the boundary - conditions of the unconfined aquifer are listed under [Constant Head](@ref) and [Boundary - conditions](@ref) of the Model parameters section. + `sbm` model type. Parameters that are part of this component are described in the [Lateral + subsurface flow](@ref params_ssf) section of Model parameters. Input parameters for this + component are derived from the SBM vertical concept and the land slope. One external + parameter [`ksathorfrac`](@ref params_ssf) is used to calculate the horizontal hydraulic + conductivity at the soil surface `kh_0`. +2. Groundwater flow (see section [Groundwater flow component](@ref lateral_gwf)) as part of the + `sbm_gwf` model type. For the unconfined aquifer the input parameters are described in the + section [Unconfined aquifer](@ref) of Model parameters. The bottom (`bottom`) of the + groundwater layer is derived from from the `soilthickness` [mm] parameter of `SBM` and the + provided surface elevation `altitude` [m] as part of the static input. The `area` parameter + is derived from the model grid. Parameters that are part of the boundary conditions of the + unconfined aquifer are listed under [Constant Head](@ref) and [Boundary conditions](@ref) of + the Model parameters section. Most hydrological model configurations make use of the kinematic wave surface routing (river -flow, overland flow or both) and input data required for the river and overland flow -components is described in [Surface flow](@ref). There is also the option to use the local -inertial model as part of the wflow\_sbm models (model types `sbm` and `sbm_gwf`): +flow, overland flow or both) and input data required for the river and overland flow components +is described in [Surface flow](@ref). There is also the option to use the local inertial model +as part of the wflow\_sbm models (model types `sbm` and `sbm_gwf`): + for river flow, see also the [Local inertial river and floodplain](@ref config_sbm_gwf_lie_river) model. + for 1D river flow and 2D overland flow combined, see also the [Local inertial river (1D) @@ -105,25 +103,25 @@ Reservoirs or lakes can be part of the kinematic wave or local inertial model fo and input parameters are described in [Reservoirs](@ref reservoir_params) and [Lakes](@ref lake_params). -The [wflow\_sediment](@ref config_sediment) model configuration consists of the vertical -[Soil Erosion](@ref) concept and the input parameters for this concept are described in the -[Sediment](@ref params_sediment) section of the Model parameters. The parameters of the -lateral [Sediment Flux in overland flow](@ref) concept are described in the [Overland -flow](@ref) section of the Model parameters. Parameters of this component are not directly -set by data from static input. The input parameters of the lateral concept [River Sediment -Model](@ref) are listed in [River flow](@ref) of the Model parameters section. +The [wflow\_sediment](@ref config_sediment) model configuration consists of the vertical [Soil +Erosion](@ref) concept and the input parameters for this concept are described in the +[Sediment](@ref params_sediment) section of the Model parameters. The parameters of the lateral +[Sediment Flux in overland flow](@ref) concept are described in the [Overland flow](@ref) +section of the Model parameters. Parameters of this component are not directly set by data from +static input. The input parameters of the lateral concept [River Sediment Model](@ref) are +listed in [River flow](@ref) of the Model parameters section. -The Model parameters section lists all the parameters per Model component and these Tables -can also be used to check which parameters can be part of the output, see also [Output -netCDF section](@ref) and [Output CSV section](@ref). +The Model parameters section lists all the parameters per Model component and these Tables can +also be used to check which parameters can be part of the output, see also [Output netCDF +section](@ref) and [Output CSV section](@ref). Example models can be found in the [Example models section](@ref sample_data). ## References -+ Yamazaki, D., Ikeshima, D., Sosa, J., Bates, P. D., Allen, G. H. and Pavelsky, T. M.: - MERIT Hydro: A high‐resolution global hydrography map based on latest topography datasets, - Water Resour. Res., 2019WR024873, doi:10.1029/2019WR024873, 2019. -+ Eilander, D., van Verseveld, W., Yamazaki, D., Weerts, A., Winsemius, H. C., and Ward, P. - J.: A hydrography upscaling method for scale-invariant parametrization of distributed ++ Yamazaki, D., Ikeshima, D., Sosa, J., Bates, P. D., Allen, G. H. and Pavelsky, T. M.: MERIT + Hydro: A high‐resolution global hydrography map based on latest topography datasets, Water + Resour. Res., 2019WR024873, doi:10.1029/2019WR024873, 2019. ++ Eilander, D., van Verseveld, W., Yamazaki, D., Weerts, A., Winsemius, H. C., and Ward, P. J.: + A hydrography upscaling method for scale-invariant parametrization of distributed hydrological models, Hydrol. Earth Syst. Sci., 25, 5287–5313, , 2021. diff --git a/docs/getting_started/download_example_models.qmd b/docs/getting_started/download_example_models.qmd index a09dc33ee..0c48853a5 100644 --- a/docs/getting_started/download_example_models.qmd +++ b/docs/getting_started/download_example_models.qmd @@ -3,8 +3,8 @@ title: Download example models --- For each wflow Model a test model is available that can help to understand the data -requirements and the usage of each Model. The TOML configuration file per available model -are listed in the Table below: +requirements and the usage of each Model. The TOML configuration file per available model are +listed in the Table below: | model | TOML configuration file | |:--------------- | ------------------| @@ -12,9 +12,9 @@ are listed in the Table below: | wflow\_sbm + groundwater flow | [sbm\_gwf\_config.toml](https://raw.githubusercontent.com/Deltares/Wflow.jl/master/test/sbm_gwf_config.toml) | | wflow_sediment | [sediment_config.toml](https://raw.githubusercontent.com/Deltares/Wflow.jl/master/test/sediment_config.toml) | -The associated Model files (input static, forcing and state files) can easily be downloaded -and for this we share the following Julia code (per Model) that downloads the required files -to your current working directory. For running these test model see also [Usage](@ref run_wflow) +The associated Model files (input static, forcing and state files) can easily be downloaded and +for this we share the following Julia code (per Model) that downloads the required files to +your current working directory. For running these test model see also [Usage](@ref run_wflow) and [Command Line Interface](@ref cli). ## wflow\_sbm + kinematic wave diff --git a/docs/getting_started/install.qmd b/docs/getting_started/install.qmd index dd5e1a80a..9406d16a6 100644 --- a/docs/getting_started/install.qmd +++ b/docs/getting_started/install.qmd @@ -1,16 +1,6 @@ --- title: Installing wflow --- -First download and install the [current stable release of -Julia](https://julialang.org/downloads/#current_stable_release). Please see [platform -specific instructions](https://julialang.org/downloads/platform/) for further installation -instructions and if you have trouble installing Julia. - -If you are new to Julia, it might be a good idea to check out the [Getting Started section -of the Julia Manual](https://docs.julialang.org/en/v1/manual/getting-started/). You can also -find additional learning resources at -[julialang.org/learning](https://julialang.org/learning/). - Wflow can be used in two different ways, depending on the required use of the code: @@ -29,14 +19,22 @@ Wflow is a [Julia](https://julialang.org/) package that can be installed in seve Below, we show how to install wflow from Julia's package repository and how to install the latest version from GitHub. +First, download and install the [current stable release of +Julia](https://julialang.org/downloads/#current_stable_release). If you have any issues +installing Julia, please see [platform specific +instructions](https://julialang.org/downloads/platform/) for further instructions. + +If you are new to Julia, it might be a good idea to check out the [Getting Started section of +the Julia Manual](https://docs.julialang.org/en/v1/manual/getting-started/). You can also find +additional learning resources at [julialang.org/learning](https://julialang.org/learning/). + ### Install from Julia's package repository To access Julia's package manager, press `]` in the Julia REPL. To get back to the Julia REPL, press backspace or ^C. ::: {.callout-tip} -If you have not used Julia in a while, it can be a good idea to run `up` to update your -packages. +If you haven't used Julia in a while, it's a good idea to run `up` to update your packages. ```julia-repl pkg> up ``` @@ -60,20 +58,20 @@ pkg> add Wflow#master This command tracks the `master` branch and updates to the latest commit on that branch when you run `update`, or simply `up`, in the Pkg REPL. The `add` installs wflow in your home -directory under `.julia/packages/Wflow`. Note that packages installed under `packages` by -`add` should not be changed in the directory, as the change could disrupt Pkg's automatic -dependency handling. +directory under `.julia/packages/Wflow`. Note that packages installed under `packages` by `add` +should not be changed in the directory, as the change could disrupt Pkg's automatic dependency +handling. -If you want to modify any files in the repository, you need to do a development install. -This can be done using: +If you want to modify any files in the repository, you need to do a development install. This +can be done using: ```julia-repl pkg> dev Wflow ``` -This will clone the git repository, place it under your home directory in -`.julia/dev/Wflow`, and add the wflow package to your project environment. To receive -updates, you'll need to pull the latest changes manually using `git pull`. +This will clone the git repository, place it under your home directory in `.julia/dev/Wflow`, +and add the wflow package to your project environment. To receive updates, you'll need to pull +the latest changes manually using `git pull`. ### Check installation of wflow @@ -84,35 +82,35 @@ julia> using Wflow ``` The first time you do this, it may take longer as any new or changed packages need to be -precompiled to enable faster loading on subsequent uses. No error messages should appear, -which indicates that you have successfully installed wflow. +precompiled to enable faster loading on subsequent uses. No error messages should appear, which +indicates that you have successfully installed wflow. Before concluding this section, we recommend a few tools that can make using and developing Julia code easier. ::: {.callout-tip} -There is a section on editors and IDEs for Julia on , scroll -down to see it. We use and recommend Microsoft's free and open source [Visual Studio +There is a section on editors and IDEs for Julia on , scroll down to +see it. We use and recommend Microsoft's free and open source [Visual Studio Code](https://code.visualstudio.com/). When combined with the [Julia -extension](https://www.julia-vscode.org/) it provides a powerful and interactive -development experience. +extension](https://www.julia-vscode.org/) it provides a powerful and interactive development +experience. ::: ::: {.callout-tip} -When planning to make changes to the code of wflow, we recommend installing the `Revise.jl` -package. This package allows you to modify code and use the changes without restarting -Julia. Install it with `add Revise` from the Pkg REPL. Then create a file called -`.julia/config/startup.jl`, and put `using Revise` there. This will load Revise every -time you start a Julia session. +If you plan to modify the code of wflow, we recommend installing the `Revise.jl` +package. This package allows you to modify code and use the changes without restarting Julia. +Install it with `add Revise` from the Pkg REPL. Then create a file called +`.julia/config/startup.jl`, and put `using Revise` there. This will load Revise every time you +start a Julia session. ::: ## Installing the compiled executable Binaries of `wflow_cli` can be downloaded from our website [download.deltares.nl](https://download.deltares.nl/en/download/wflow/), and are currently -available for Windows. Download and install the `.msi` file. After installation, you will -see two folders in the installation directory. Only the `bin/wflow_cli` is used. The -`artifacts` folder contains binary dependencies such as netCDF. +available for Windows. Download and install the `.msi` file. After installation, you will see +two folders in the installation directory. Only the `bin/wflow_cli` is used. The `artifacts` +folder contains binary dependencies such as netCDF. ``` artifacts\ @@ -128,10 +126,10 @@ Usage: wflow_cli 'path/to/config.toml' ::: {.callout-note} -The old version of wflow, based on Python and PCRaster libraries is also available to -download from our website [download.deltares.nl](https://download.deltares.nl/en/download/wflow/). -We recommend installing the Julia version, as this documentation is written to support -this version. +The old version of wflow, based on Python and PCRaster libraries is also available to download +from our website [download.deltares.nl](https://download.deltares.nl/en/download/wflow/). We +recommend installing the Julia version, as this documentation is written to support this +version. ::: diff --git a/docs/getting_started/running_wflow.qmd b/docs/getting_started/running_wflow.qmd index 565846c99..b59194a4f 100644 --- a/docs/getting_started/running_wflow.qmd +++ b/docs/getting_started/running_wflow.qmd @@ -4,8 +4,8 @@ title: Running a simulation ## Using Julia -Once you installed Julia and Wflow.jl, a simulation can be started from the command line -as follows: +Once you installed Julia and Wflow.jl, a simulation can be started from the command line as +follows: ```bash julia -e 'using Wflow; Wflow.run()' path/to/config.toml @@ -29,20 +29,19 @@ config.endtime = DateTime("2000-01-03T00:00:00") Wflow.run(config) ``` -## [Using the command line interface](@id cli) +## Using the command line interface If you don't need the extra features of using wflow as a library, but just want to run simulations, the command line interface makes it easier to do so. It consists of a single -executable, `wflow_cli` that accepts a single argument, the path to a TOML configuration -file. +executable, `wflow_cli` that accepts a single argument, the path to a TOML configuration file. Binaries of `wflow_cli` can be downloaded from our website [download.deltares.nl](https://download.deltares.nl/en/download/wflow/), and are currently available for Windows. After installing you can see three folders in the installation directory. It is only the -`bin/wflow_cli` that is directly used. All three folders need to stay together however. -The share folder contains TOML files with more information about the build. +`bin/wflow_cli` that is directly used. All three folders need to stay together however. The +share folder contains TOML files with more information about the build. ```bash bin\wflow_cli diff --git a/docs/home/case_studies.qmd b/docs/home/case_studies.qmd index be5a3debe..1367bf3d2 100644 --- a/docs/home/case_studies.qmd +++ b/docs/home/case_studies.qmd @@ -5,78 +5,75 @@ title: "Case studies" ## Wflow models for the Meuse and Rhine Reliable hydrological models for the Rhine and the Meuse river basins are necessary for -short-term forecasting of river flows and long-term predictions for strategic water -management planning. In collaboration with Rijkswaterstaat, Deltares is developing a new -line of models for the Rhine and the Meuse basins. The models will be used for forecasting -and to estimate the impact of climate change on water resources and extreme streamflow. In -the model development, we aim to improve hydrological predictions by including relevant -processes in the model schematization. The modularity of the wflow framework is ideal for -this as we can easily evaluate the combination of different vertical and lateral model -components. For example, the local inertial routing for river and overland flow enables us -to consider retention of water in the floodplains, which is likely to improve extreme -streamflow predictions. +short-term forecasting of river flows and long-term predictions for strategic water management +planning. In collaboration with Rijkswaterstaat, Deltares is developing a new line of models +for the Rhine and the Meuse basins. The models will be used for forecasting and to estimate the +impact of climate change on water resources and extreme streamflow. In the model development, +we aim to improve hydrological predictions by including relevant processes in the model +schematization. The modularity of the wflow framework is ideal for this as we can easily +evaluate the combination of different vertical and lateral model components. For example, the +local inertial routing for river and overland flow enables us to consider retention of water in +the floodplains, which is likely to improve extreme streamflow predictions. ![](../images/case_rhine_meuse.png) ## Operational flood forecasting in Australia -In Australia, there was a need for high-resolution, fast and accurate rainfall-runoff models -to provide boundary conditions for a fast and detailed flood inundation model (SFINCS). The -domain of the flood model covers the entire North and East Coast of Australia. Although -many gauging stations are available to provide real-time information, many rivers are not -covered. For these locations, wflow\_sbm models are used to provide this real-time -information. Additionally, these models are used to provide projections for potential future -scenarios. Using the HydroMT library, all wflow\_sbm models were automatically built. The -high level of flexibility in spatial and temporal resolution, combined with the physics-based nature -of the concept, makes Wflow\_sbm particularly suitable for ungauged basins. Furthermore, the -model is detailed and computationally efficient enough for coupling with the fast -flood inundation model SFINCS. +In Australia, there was a need for high-resolution, fast and accurate rainfall-runoff models to +provide boundary conditions for a fast and detailed flood inundation model (SFINCS). The domain +of the flood model covers the entire North and East Coast of Australia. Although many gauging +stations are available to provide real-time information, many rivers are not covered. For these +locations, wflow\_sbm models are used to provide this real-time information. Additionally, +these models are used to provide projections for potential future scenarios. Using the HydroMT +library, all wflow\_sbm models were automatically built. The high level of flexibility in +spatial and temporal resolution, combined with the physics-based nature of the concept, makes +Wflow\_sbm particularly suitable for ungauged basins. Furthermore, the model is detailed and +computationally efficient enough for coupling with the fast flood inundation model SFINCS. ![](../images/case_flifs_1.png) -The results of this proof of concept are very promising. Technically, we were able to -quickly set up the wflow\_sbm models, couple them to the flood inundation models (SFINCS), and -run the models operationally under the Delft-FEWS platform. Model validation was conducted -for two basins by comparing the results of Wflow\_sbm against observations and the -results of calibrated URBS models. This validation demonstrated that the uncalibrated Wflow\_sbm -model results were already quite satisfactory, especially given the complex nature of these -basins, which include several small and large reservoirs. We could also show the potential -for further calibration by adjusting the KsatHorFrac parameter. - -Reference: De Kleermaeker, S., Leijnse, T., Morales, Y., Druery, C., Maguire, S., -1. Developing a real-time data and modelling framework for operational flood inundation - forecasting in Australia. In Hydrology & Water Resources Symposium 2022 (HWRS 2022): - The Past, the Present, the Future. Engineers Australia. - https://search.informit.org/doi/10.3316/informit.916755150845355 +The results of this proof of concept are very promising. Technically, we were able to quickly +set up the wflow\_sbm models, couple them to the flood inundation models (SFINCS), and run the +models operationally under the Delft-FEWS platform. Model validation was conducted for two +basins by comparing the results of Wflow\_sbm against observations and the results of +calibrated URBS models. This validation demonstrated that the uncalibrated Wflow\_sbm model +results were already quite satisfactory, especially given the complex nature of these basins, +which include several small and large reservoirs. We could also show the potential for further +calibration by adjusting the KsatHorFrac parameter. + +De Kleermaeker, S., Leijnse, T., Morales, Y., Druery, C., & Maguire, S. (2022). Developing a +real-time data and modelling framework for operational flood inundation forecasting in +Australia. In Hydrology & Water Resources Symposium 2022 (HWRS 2022): The Past, the Present, +the Future. Engineers Australia. +https://search.informit.org/doi/10.3316/informit.916755150845355 ![](../images/case_flifs_2.png) ## Simulating plastic transport in Thailand -For the Pollution Control Board of the Government of Thailand and the World Bank, we -supported a material flow analysis of plastics in Thailand using wflow. Plastic pollution -is a growing global issue. Plastic waste enters rivers and is transported to the ocean -where it persists and threatens the health of the ocean, seas and coasts. The initial -movement of plastic waste is in many cases triggered by runoff from (heavy) rainfall and -transported by water flow towards small streams and rivers. Therefore there is strong -relation to rainfall-runoff processes, which can be modeled using high-resolution -rainfall-runoff models. +For the Pollution Control Board of the Government of Thailand and the World Bank, we supported +a material flow analysis of plastics in Thailand using wflow. Plastic pollution is a growing +global issue. Plastic waste enters rivers and is transported to the ocean where it persists and +threatens the health of the ocean, seas and coasts. The initial movement of plastic waste is in +many cases triggered by runoff from (heavy) rainfall and transported by water flow towards +small streams and rivers. Therefore there is strong relation to rainfall-runoff processes, +which can be modeled using high-resolution rainfall-runoff models. In this study we applied the wflow\_sbm model in combination with a fate-and-transport and water quality model (DelWaq) to simulate the movement of plastics through five large river -basins and on three island and coastal zones (Krabi, Phuket, and Ko Samui; see screenshot of the -model below) in Thailand. Together with our partners Panya Consultants and HII, we were able -to identify hotspots of plastic pollution, estimate how much plastic waste would end up in the Gulf of -Thailand and recommend priority areas for reducing plastic waste reaching the sea. +basins and on three island and coastal zones (Krabi, Phuket, and Ko Samui; see screenshot of +the model below) in Thailand. Together with our partners Panya Consultants and HII, we were +able to identify hotspots of plastic pollution, estimate how much plastic waste would end up in +the Gulf of Thailand and recommend priority areas for reducing plastic waste reaching the sea. ![](../images/case_mfa_1.png) The wflow\_sbm models for the five large basins were calibrated. The presence of large dams and -reservoirs complicated calibration, but with the input for the dam operation, -the model performance for these basins could be largely improved. The figure below shows the -calibrated model results for the Chao Phraya, just upstream of Bangkok. The input from the -hydrological wflow\_sbm model was used as input for the fate and transport model to assess -the amount of plastic transported to the ocean. +reservoirs complicated calibration, but with the input for the dam operation, the model +performance for these basins could be largely improved. The figure below shows the calibrated +model results for the Chao Phraya, just upstream of Bangkok. The input from the hydrological +wflow\_sbm model was used as input for the fate and transport model to assess the amount of +plastic transported to the ocean. ![](../images/case_mfa_3.png) diff --git a/docs/home/publications.qmd b/docs/home/publications.qmd index 1e1cf6a26..a05714cac 100644 --- a/docs/home/publications.qmd +++ b/docs/home/publications.qmd @@ -7,34 +7,34 @@ title: "Publications" For publications, please cite the following paper introducing Wflow.jl and describing the wflow\_sbm concept, together with some case studies: -van Verseveld, W. J., Weerts, A. H., Visser, M., Buitink, J., Imhoff, R. O., Boisgontier, -H., Bouaziz, L., Eilander, D., Hegnauer, M., ten Velden, C., and Russell, B., 2024. -Wflow_sbm v0.7.3, a spatially distributed hydrological model: from global data to local -applications. Geosci. Model Dev., 17, 3199–3234. . +van Verseveld, W. J., Weerts, A. H., Visser, M., Buitink, J., Imhoff, R. O., Boisgontier, H., +Bouaziz, L., Eilander, D., Hegnauer, M., ten Velden, C., and Russell, B., 2024. Wflow_sbm +v0.7.3, a spatially distributed hydrological model: from global data to local applications. +Geosci. Model Dev., 17, 3199–3234. . To cite a specific software version please use the DOI provided in the Zenodo badge [![DOI](https://zenodo.org/badge/246787232.svg)](https://zenodo.org/badge/latestdoi/246787232), -that points to the latest release. The DOIs of previous versions are also available at -Zenodo. If you use a snapshot of the development (without a DOI) please cite as follows: +that points to the latest release. The DOIs of previous versions are also available at Zenodo. +If you use a snapshot of the development (without a DOI) please cite as follows: -van Verseveld, Willem, Visser, Martijn, Buitink, Joost, Bouaziz, Laurène, Boisgontier, -Hélène, Bootsma, Huite, Weerts, Albrecht, Baptista, Carlos Fernando, Pronk, Maarten, -Eilander, Dirk, Hartgring, Sebastian, Dalmijn, Brendan, Hofer, Julian, Hegnauer, Mark, & -Mendoza, Raul, (YEAR). Deltares/Wflow.jl: unstable-master. -, obtained: DATE\_OF\_DOWNLOAD. +van Verseveld, Willem, Visser, Martijn, Buitink, Joost, Bouaziz, Laurène, Boisgontier, Hélène, +Bootsma, Huite, Weerts, Albrecht, Baptista, Carlos Fernando, Pronk, Maarten, Eilander, Dirk, +Hartgring, Sebastian, Dalmijn, Brendan, Hofer, Julian, Hegnauer, Mark, & Mendoza, Raul, (YEAR). +Deltares/Wflow.jl: unstable-master. , obtained: +DATE\_OF\_DOWNLOAD. ## Publications using wflow ### Peer reviewed journal papers -Aerts, J. P. M., Hut, R. W., van de Giesen, N. C., Drost, N., van Verseveld, W. J., Weerts, -A. H., and Hazenberg, P., 2022. Large-sample assessment of varying spatial resolution on the +Aerts, J. P. M., Hut, R. W., van de Giesen, N. C., Drost, N., van Verseveld, W. J., Weerts, A. +H., and Hazenberg, P., 2022. Large-sample assessment of varying spatial resolution on the streamflow estimates of the wflow_sbm hydrological model. Hydrol. Earth Syst. Sci., 26, 4407–4430. . -de Boer-Euser, T., Bouaziz, L., De Niel, J., Brauer, C., Dewals, B., Drogue, G., Fenicia, -F., Grelier, B., Nossent, J., Pereira, F., Savenije, H., Thirel, G., Willems, P., 2017. -Looking beyond general metrics for model comparison – lessons from an international model +de Boer-Euser, T., Bouaziz, L., De Niel, J., Brauer, C., Dewals, B., Drogue, G., Fenicia, F., +Grelier, B., Nossent, J., Pereira, F., Savenije, H., Thirel, G., Willems, P., 2017. Looking +beyond general metrics for model comparison – lessons from an international model intercomparison study. Hydrol. Earth Syst. Sci. 21, 423–440. . @@ -44,44 +44,43 @@ Pereira, F., Sprokkereef, E., Stam, J., Weerts, A.H., Willems, P., Savenije, H.H Hrachowitz, M., 2021. Behind the scenes of streamflow model performance. Hydrol. Earth Syst. Sci., 25, 1069–1095. . -Bouaziz, L. J. E., Aalbers, E. E., Weerts, A. H., Hegnauer, M., Buiteveld, H., Lammersen, -R., Stam, J., Sprokkereef, E., Savenije, H. H. G., and Hrachowitz, M., 2022. Ecosystem -adaptation to climate change: the sensitivity of hydrological predictions to time-dynamic -model parameters, Hydrol. Earth Syst. Sci., 26, 1295–1318. +Bouaziz, L. J. E., Aalbers, E. E., Weerts, A. H., Hegnauer, M., Buiteveld, H., Lammersen, R., +Stam, J., Sprokkereef, E., Savenije, H. H. G., and Hrachowitz, M., 2022. Ecosystem adaptation +to climate change: the sensitivity of hydrological predictions to time-dynamic model +parameters, Hydrol. Earth Syst. Sci., 26, 1295–1318. . -Casson, D. R., Werner, M., Weerts, A., and Solomatine, D., 2018. Global re-analysis datasets -to improve hydrological assessment and snow water equivalent estimation in a sub-Arctic -watershed. Hydrol. Earth Syst. Sci., 22, 4685–4697. -. +Casson, D. R., Werner, M., Weerts, A., and Solomatine, D., 2018. Global re-analysis datasets to +improve hydrological assessment and snow water equivalent estimation in a sub-Arctic watershed. +Hydrol. Earth Syst. Sci., 22, 4685–4697. . -Droppers, B., Rakovec, O., Avila, L., Azimi, S., Cortés-Torres N., De León Pérez, D., -Imhoff, R., Francés, F., Kollet, S., Rigon, R., Weerts, A. & Samaniego, L, 2024. Multi-model +Droppers, B., Rakovec, O., Avila, L., Azimi, S., Cortés-Torres N., De León Pérez, D., Imhoff, +R., Francés, F., Kollet, S., Rigon, R., Weerts, A. & Samaniego, L, 2024. Multi-model hydrological reference dataset over continental Europe and an African basin. Sci Data, 11, 1009\. . Emerton, R.E., Stephens, E.M., Pappenberger, F., Pagano, T.C., Weerts, A.H., Wood, A.W., -Salamon, P., Brown, J.D., Hjerdt, N., Donnelly, C., Baugh, C.A., Cloke, H.L., 2016. -Continental and global scale flood forecasting systems. WIREs Water 3, 391–418. +Salamon, P., Brown, J.D., Hjerdt, N., Donnelly, C., Baugh, C.A., Cloke, H.L., 2016. Continental +and global scale flood forecasting systems. WIREs Water 3, 391–418. . -Gebremicael, T.G., Mohamed, Y.A., Van der Zaag, P., 2019. Attributing the hydrological -impact of different land use types and their long-term dynamics through combining -parsimonious hydrological modelling, alteration analysis and PLSR analysis. Science of The -Total Environment, 660, 1155-1167, . +Gebremicael, T.G., Mohamed, Y.A., Van der Zaag, P., 2019. Attributing the hydrological impact +of different land use types and their long-term dynamics through combining parsimonious +hydrological modelling, alteration analysis and PLSR analysis. Science of The Total +Environment, 660, 1155-1167, . Giardino, A., Schrijvershof, R., Nederhoff, C.M., de Vroeg, H., Brière, C., Tonnon, P.-K., -Caires, S., Walstra, D.J., Sosa, J., van Verseveld, W., Schellekens, J., Sloff, C.J., 2018. -A quantitative assessment of human interventions and climate change on the West African -sediment budget, Ocean & Coastal Management, 156, 249-265. +Caires, S., Walstra, D.J., Sosa, J., van Verseveld, W., Schellekens, J., Sloff, C.J., 2018. A +quantitative assessment of human interventions and climate change on the West African sediment +budget, Ocean & Coastal Management, 156, 249-265. . Hally, A., Caumont, O., Garrote, L., Richard, E., Weerts, A., Delogu, F., Fiori, E., Rebora, N., Parodi, A., Mihalović, A., Ivković, M., Dekić, L., van Verseveld, W., Nuissier, O., -Ducrocq, V., D’Agostino, D., Galizia, A., Danovaro, E., Clematis, A., 2015. -Hydrometeorological multi-model ensemble simulations of the 4 November 2011 flash flood -event in Genoa, Italy, in the framework of the DRIHM project. Nat. Hazards Earth Syst. Sci. -15, 537–555. . +Ducrocq, V., D’Agostino, D., Galizia, A., Danovaro, E., Clematis, A., 2015. Hydrometeorological +multi-model ensemble simulations of the 4 November 2011 flash flood event in Genoa, Italy, in +the framework of the DRIHM project. Nat. Hazards Earth Syst. Sci. 15, 537–555. +. Hassaballah, K., Mohamed, Y., Uhlenbrook, S., and Biro, K., 2017. Analysis of streamflow response to land use and land cover changes using satellite data and hydrological modelling: @@ -89,18 +88,18 @@ case study of Dinder and Rahad tributaries of the Blue Nile (Ethiopia–Sudan), Syst. Sci., 21, 5217–5242. . Imhoff, R.O., Buitink, J., van Verseveld, W.J., Weerts, A.H., 2024. A fast high resolution -distributed hydrological model for forecasting, climate scenarios and digital twin -applications using wflow_sbm. Environmental Modelling & Software, 179, 106099. +distributed hydrological model for forecasting, climate scenarios and digital twin applications +using wflow_sbm. Environmental Modelling & Software, 179, 106099. -Imhoff, R.O, van Verseveld, W.J., van Osnabrugge, B., Weerts, A.H., 2020. Scaling -Point-Scale (Pedo)transfer Functions to Seamless Large-Domain Parameter Estimates for -High-Resolution Distributed Hydrologic Modeling: An Example for the Rhine River. Water -Resources Research, 56, e2019WR026807. . +Imhoff, R.O, van Verseveld, W.J., van Osnabrugge, B., Weerts, A.H., 2020. Scaling Point-Scale +(Pedo)transfer Functions to Seamless Large-Domain Parameter Estimates for High-Resolution +Distributed Hydrologic Modeling: An Example for the Rhine River. Water Resources Research, 56, +e2019WR026807. . -Imhoff, R.O., van Verseveld, W., van Osnabrugge, B., Weerts, A.H., 2020. Ruimtelijk -schaalbare hydrologische modelparameters uit open-source omgevingsdata: een voorbeeld voor -de Rijn. Stromingen: vakblad voor hydrologen, 26(3), 19-36 . +Imhoff, R.O., van Verseveld, W., van Osnabrugge, B., Weerts, A.H., 2020. Ruimtelijk schaalbare +hydrologische modelparameters uit open-source omgevingsdata: een voorbeeld voor de Rijn. +Stromingen: vakblad voor hydrologen, 26(3), 19-36 . Jeuken, A., Bouaziz, L., Corzo, G., Alfonso, L., 2016. Analyzing Needs for Climate Change Adaptation in the Magdalena River Basin in Colombia, in: Filho, W.L., Musa, H., Cavan, G., @@ -110,13 +109,13 @@ Change Management. Springer International Publishing, pp. 329–344 López López, P., Wanders, N., Schellekens, J., Renzullo, L.J., Sutanudjaja, E.H., Bierkens, M.F.P., 2016. Improved large-scale hydrological modelling through the assimilation of -streamflow and downscaled satellite soil moisture observations. Hydrol. Earth Syst. Sci., -20, 3059–3076. . +streamflow and downscaled satellite soil moisture observations. Hydrol. Earth Syst. Sci., 20, +3059–3076. . -Pranoto, B., Soekarno, H., Hartulistiyoso, E., Nur Aidi, M., Sutrisno, D., Pohan, D., -Radhika, Sutejo, B., Heru Kuncoro, A., Nahib, I., 2024. Integrating Flood Early Warning -System (FEWS) for Optimizing Small Hydropower Sites: A West Java Case Study. EVERGREEN Joint -Journal of Novel Carbon Resource Sciences & Green Asia Strategy, 11, 3, 2691-2699. +Pranoto, B., Soekarno, H., Hartulistiyoso, E., Nur Aidi, M., Sutrisno, D., Pohan, D., Radhika, +Sutejo, B., Heru Kuncoro, A., Nahib, I., 2024. Integrating Flood Early Warning System (FEWS) +for Optimizing Small Hydropower Sites: A West Java Case Study. EVERGREEN Joint Journal of Novel +Carbon Resource Sciences & Green Asia Strategy, 11, 3, 2691-2699. Rakovec, O., Weerts, A.H., Sumihar, J., Uijlenhoet, R., 2015. Operational aspects of @@ -124,49 +123,49 @@ asynchronous filtering for flood forecasting. Hydrol. Earth Syst. Sci., 19, 2911 . Ratri, D.N., A.H. Weerts, R. Muharsyah, K. Whan, A. Klein Tank, E. Aldrian, M. H. Hariadi, -Calibration of ECMWF SEAS5 based streamflow forecast in Seasonal hydrological forecasting -for Citarum river basin, West Java, Indonesia, Journal of Hydrology: Regional -Studies,45,101305, . +Calibration of ECMWF SEAS5 based streamflow forecast in Seasonal hydrological forecasting for +Citarum river basin, West Java, Indonesia, Journal of Hydrology: Regional Studies,45,101305, +. Rusli, S.R., A.H. Weerts, A. Taufiq, V. Bense, 2021. Estimating water balance components and -their uncertainty bounds in highly groundwater-dependent and data-scarce area: An example -for the Upper Citarum basin, J. Hydrol. Regional Studies, +their uncertainty bounds in highly groundwater-dependent and data-scarce area: An example for +the Upper Citarum basin, J. Hydrol. Regional Studies, . Rusli, S.R., V.F. Bense, A. Taufiq, A.H. Weerts,2023. Quantifying basin-scale changes in -groundwater storage using GRACE and one-way coupled hydrological and groundwater flow model -in the data-scarce Bandung groundwater Basin, Indonesia, Groundwater for Sustainable +groundwater storage using GRACE and one-way coupled hydrological and groundwater flow model in +the data-scarce Bandung groundwater Basin, Indonesia, Groundwater for Sustainable Development,22, 100953, . -Rusli, S.R., A.H. Weerts, S.M.T. Mustafa, D.E. Irawan, A. Taufiq, V.F. Bense, 2023. -Quantifying aquifer interaction using numerical groundwater flow model evaluated by -environmental water tracer data: Application to the data-scarce area of the Bandung -groundwater basin, West Java, Indonesia, Journal of Hydrology: Regional Studies, 50, +Rusli, S.R., A.H. Weerts, S.M.T. Mustafa, D.E. Irawan, A. Taufiq, V.F. Bense, 2023. Quantifying +aquifer interaction using numerical groundwater flow model evaluated by environmental water +tracer data: Application to the data-scarce area of the Bandung groundwater basin, West Java, +Indonesia, Journal of Hydrology: Regional Studies, 50, . -Schaller, N., Sillmann, J., Müller, M., Haarsma, R., Hazeleger, W., Hegdahl, T.J., Kelder, -T., van den Oord, G., Weerts, A., Whan, K., 2020. The role of spatial and temporal model -resolution in a flood event storyline approach in western Norway. Weather and Climate -Extremes, doi: . +Schaller, N., Sillmann, J., Müller, M., Haarsma, R., Hazeleger, W., Hegdahl, T.J., Kelder, T., +van den Oord, G., Weerts, A., Whan, K., 2020. The role of spatial and temporal model resolution +in a flood event storyline approach in western Norway. Weather and Climate Extremes, doi: +. Seizarwati, W. and M. Syahidah, 2021. Rainfall-Runoff Simulation for Water Availability -Estimation in Small Island Using Distributed Hydrological Model wflow. IOP Conf. Ser.: -Earth Environ. Sci., 930,012050, doi:10.1088/1755-1315/930/1/012050. +Estimation in Small Island Using Distributed Hydrological Model wflow. IOP Conf. Ser.: Earth +Environ. Sci., 930,012050, doi:10.1088/1755-1315/930/1/012050. -Sperna Weiland, F.C., R.D. Visser, P. Greve, B. Bisselink, L. Brunner and A.H. Weerts, -2021\. Estimating Regionalized Hydrological Impacts of Climate Change Over Europe by -Performance-Based Weighting of CORDEX Projections, Frontiers of Water, +Sperna Weiland, F.C., R.D. Visser, P. Greve, B. Bisselink, L. Brunner and A.H. Weerts, 2021\. +Estimating Regionalized Hydrological Impacts of Climate Change Over Europe by Performance-Based +Weighting of CORDEX Projections, Frontiers of Water, . Tangdamrongsub, N., Steele-Dunne, S.C., Gunter, B.C., Ditmar, P.G., Weerts, A.H., 2015. Data -assimilation of GRACE terrestrial water storage estimates into a regional hydrological model -of the Rhine River basin. Hydrol. Earth Syst. Sci. 19, 2079–2100. +assimilation of GRACE terrestrial water storage estimates into a regional hydrological model of +the Rhine River basin. Hydrol. Earth Syst. Sci. 19, 2079–2100. . -van der Laan, E., P. Hazenberg, A.H. Weerts, 2024. Simulation of long-term storage dynamics -of headwater reservoirs across the globe using public cloud computing infrastructure. -Science of The Total Environment, 172678, . +van der Laan, E., P. Hazenberg, A.H. Weerts, 2024. Simulation of long-term storage dynamics of +headwater reservoirs across the globe using public cloud computing infrastructure. Science of +The Total Environment, 172678, . van Osnabrugge, B., Weerts, A.H., Uijlenhoet, R., 2017. genRE: A method to extend gridded precipitation climatology data sets in near real-time for hydrological forecasting purposes. @@ -178,35 +177,34 @@ forecasts to 10-day streamflow forecast skill for the Rhine River, Hydrol. Earth van der Vat, M., Boderie, P., Bons, K.A., Hegnauer, M., Hendriksen, G., van Oorschot, M., Ottow, B., Roelofsen, F., Sankhua, R.N., Sinha, S.K., Warren, A., Young, W., 2019. -Participatory Modelling of Surface and Groundwater to Support Strategic Planning in the -Ganga Basin in India. Water, 11, 2443. . +Participatory Modelling of Surface and Groundwater to Support Strategic Planning in the Ganga +Basin in India. Water, 11, 2443. . -Wannasin, C., Brauer, C. C., Uijlenhoet, R., van Verseveld, W. J., Weerts, A. H., 2021. -Daily flow simulation in Thailand Part I: Testing a distributed hydrological model with -seamless parameter maps based on global data. Journal of Hydrology: Regional Studies, 34, -1-19. . +Wannasin, C., Brauer, C. C., Uijlenhoet, R., van Verseveld, W. J., Weerts, A. H., 2021. Daily +flow simulation in Thailand Part I: Testing a distributed hydrological model with seamless +parameter maps based on global data. Journal of Hydrology: Regional Studies, 34, 1-19. +. -Wannasin, C., Brauer, C. C., Uijlenhoet, R., van Verseveld, W. J., Weerts, A. H., 2021. -Daily flow simulation in Thailand Part II: Unraveling effects of reservoir operation. -Journal of Hydrology: RegionalStudies, 34, 1-17. -. +Wannasin, C., Brauer, C. C., Uijlenhoet, R., van Verseveld, W. J., Weerts, A. H., 2021. Daily +flow simulation in Thailand Part II: Unraveling effects of reservoir operation. Journal of +Hydrology: RegionalStudies, 34, 1-17. . Wang, X., Zhang, J., Babovic, V., 2016. Improving real-time forecasting of water quality -indicators with combination of process-based models and data assimilation technique. -Ecological Indicators 66, 428–439. . +indicators with combination of process-based models and data assimilation technique. Ecological +Indicators 66, 428–439. . ### PhD, MSc, BSc Theses & Internship reports -Abdelnour, A., 2022. Bias Correction of Climate Simulations to Assess Climate Change Impacts -on Low Flows in the Rhine River, MSc thesis, Delft Universitry of technology, +Abdelnour, A., 2022. Bias Correction of Climate Simulations to Assess Climate Change Impacts on +Low Flows in the Rhine River, MSc thesis, Delft Universitry of technology, . -Ali, M.A., 2023. Machine learning for predicting spatially variable lateral hydraulic conductivity: -a step towards efficient hydrological model calibration and global applicability, Intersnhip -report, Deltares. +Ali, M.A., 2023. Machine learning for predicting spatially variable lateral hydraulic +conductivity: a step towards efficient hydrological model calibration and global applicability, +Intersnhip report, Deltares. -Arnal, L., 2014. An intercomparison of flood forecasting models for the Meuse River basin, -MSc Thesis, Vrije Universiteit, Amsterdam, . +Arnal, L., 2014. An intercomparison of flood forecasting models for the Meuse River basin, MSc +Thesis, Vrije Universiteit, Amsterdam, . Alkemade, G., 2019. Routing and calibration of distributed hydrological models, MSc Thesis, Vrije Universiteit, Amsterdam, Faculty of Science, Hydrology. @@ -214,29 +212,28 @@ Vrije Universiteit, Amsterdam, Faculty of Science, Hydrology. Azadeh Karami Fard, 2015. Modeling runoff of an Ethiopian catchment with WFLOW, MSc. Thesis, Vrije Universiteit, Amsterdam. -Benschop, J., 2023. The application of hybrid lateral routing in hydrological simulations -of the Rhine catchment, MSc Thesis, Hydrology and Environmental Hydraulics Group, Wageningen +Benschop, J., 2023. The application of hybrid lateral routing in hydrological simulations of +the Rhine catchment, MSc Thesis, Hydrology and Environmental Hydraulics Group, Wageningen University. Beusen, B., 2021. The effect of rooting depth on discharge and evapotranspiration in (semi-)arid areas, MSc Thesis, Hydrology and Quantitative Water Management Group, Wageningen University. -Beusen, B., 2021. Plastic transport and the effect of climate change in the Rhine, -Internship report, Deltares. +Beusen, B., 2021. Plastic transport and the effect of climate change in the Rhine, Internship +report, Deltares. -Bouaziz, L. J. E., 2021. Internal processes in hydrological models: A glance at the Meuse -basin from space. Delft University of Technology, Delft, the Netherlands, Doctoral -dissertation, . +Bouaziz, L. J. E., 2021. Internal processes in hydrological models: A glance at the Meuse basin +from space. Delft University of Technology, Delft, the Netherlands, Doctoral dissertation, +. Hartgring, S., 2023. On Forecasting the Rur River: Using hindcasts and forecasts of the 2021 flood event to improve understanding of flood forecasting in the Rur catchment. Delft University of Technology, Delft, the Netherlands, MSc thesis, . -Jaime, D.E.V, 2021. Ensemble hydrological forecasts to derive extreme return -periods: Case Study of the Overijsselse Vecht River using the wflow_sbm model, MSc thesis, -Unesco IHE, Delft. +Jaime, D.E.V, 2021. Ensemble hydrological forecasts to derive extreme return periods: Case +Study of the Overijsselse Vecht River using the wflow_sbm model, MSc thesis, Unesco IHE, Delft. López López, P., 2018. Application of Global Hydrological Datasets for River Basin Modelling Utrecht University, Utrecht, the Netherlands, pp. 1-214, Doctoral dissertation, @@ -246,12 +243,12 @@ Maat, W.H., 2015. Simulating discharges and forecasting floods using a conceptua rainfall-runoff model for the Bolivian Mamoré basin, MSc Thesis, University of Twente, Enschede. . -van Osnabrugge, B., 2020. Interpolate, simulate, assimilate: operational aspects of -improving hydrological forecasts in the Rhine basin. Wageningen University, Doctoral -dissertation, . +van Osnabrugge, B., 2020. Interpolate, simulate, assimilate: operational aspects of improving +hydrological forecasts in the Rhine basin. Wageningen University, Doctoral dissertation, +. -Rohrmueller, I., 2019. BENCHMARKING THE NEW WFLOW DISTRIBUTED HYDROLOGICAL MODEL, MSc -Thesis, School of Engineering - Newcastle University. +Rohrmueller, I., 2019. BENCHMARKING THE NEW WFLOW DISTRIBUTED HYDROLOGICAL MODEL, MSc Thesis, +School of Engineering - Newcastle University. Rusli, S.R., 2024. Deepening the quantitative understanding of groundwater systems in data-scarce areas: application in the Bandung Groundwater Basin, Indonesia. Wageningen @@ -262,12 +259,11 @@ assimilation on the performance of hydrological forecasting in the Karasu catchm MSc thesis, Wageningen University. van der Gaast, R.H., 2024. Evaluating the transferability of data-driven pedo-transfer -functions for the wflow\_sbm parameter KsatHorFrac in central and Western Europe. -Universiteit Twente, Enschede, The Netherlands, -. +functions for the wflow\_sbm parameter KsatHorFrac in central and Western Europe. Universiteit +Twente, Enschede, The Netherlands, . -Verbrugge, M., 2019. Reservoir Operation Optimization, a case study in the Chao Phraya -Basin, BSc thesis, Hydrology and Quantitative Water Management Group, Wageningen University. +Verbrugge, M., 2019. Reservoir Operation Optimization, a case study in the Chao Phraya Basin, +BSc thesis, Hydrology and Quantitative Water Management Group, Wageningen University. Verbrugge, M., 2023. Bias-correcting meteorological forcing to improve seasonal discharge forecasting of the Rhine, Internship report, Deltares. @@ -280,8 +276,8 @@ Visser, B., 2020. Impact of climate change on local water resources of European Intersnhip report, Deltares. Wannasin, C., 2023. Modelling and forecasting daily streamflow with reservoir operation in the -upper Chao Phraya River basin, Thailand. -Wageningen University, Doctoral dissertation, . +upper Chao Phraya River basin, Thailand. Wageningen University, Doctoral dissertation, +. ### Reports diff --git a/docs/index.qmd b/docs/index.qmd index 2b3476bd2..9411868bb 100644 --- a/docs/index.qmd +++ b/docs/index.qmd @@ -17,35 +17,33 @@ format: # About wflow -Wflow is Deltares' solution for modelling hydrological processes, allowing users to account -for precipitation, interception, snow accumulation and melt, evapotranspiration, soil water, +Wflow is Deltares' solution for modelling hydrological processes, allowing users to account for +precipitation, interception, snow accumulation and melt, evapotranspiration, soil water, surface water, groundwater recharge, and water demand and allocation in a fully distributed environment. Successfully applied worldwide for analyzing flood hazards, drought, climate change impacts and land use changes, wflow is growing to be a leader in hydrology solutions. Wflow is conceived as a framework, within which multiple distributed model concepts are -available, which maximizes the use of open earth observation data, making it the -hydrological model of choice for data scarce environments. Based on gridded topography, -soil, land use and climate data, wflow calculates all hydrological fluxes at any given grid -cell in the model at a given time step. +available, which maximizes the use of open earth observation data, making it the hydrological +model of choice for data scarce environments. Based on gridded topography, soil, land use and +climate data, wflow calculates all hydrological fluxes at any given grid cell in the model at a +given time step. -Wflow was born out of the creation of Deltares in 2008, when a strategic review identified -the need for a distributed hydrological model to allow the simulation of flows at the -catchment scale. With the intention being to encourage greater scientific collaboration. -For this reason: +Wflow was born out of the creation of Deltares in 2008, when a strategic review identified the +need for a distributed hydrological model to allow the simulation of flows at the catchment +scale. With the intention being to encourage greater scientific collaboration. For this reason: * Wflow is free and open source software. * Wflow is easily coupled with other models and software applications. * Contribution to the wflow code development is encouraged. -From 2021 the [wflow code](https://github.com/Deltares/Wflow.jl) is distributed under the -[MIT License](https://github.com/Deltares/Wflow.jl/blob/master/LICENSE). Wflow is also -available as a [compiled executable](https://download.deltares.nl/en/download/wflow/) under -the Deltares terms and conditions. The wflow computational engine is built in the -[Julia](https://julialang.org/) language, a high-performance computing language. -Wflow does not include a graphical user interface and is designed for maximum user -flexibility. Prior to 2021, wflow was developed in Python on top of the PCRaster Python -extension. The Python version is [still available](https://github.com/openstreams/wflow), -but not actively developed. +From 2021 the [wflow code](https://github.com/Deltares/Wflow.jl) is distributed under the [MIT +License](https://github.com/Deltares/Wflow.jl/blob/master/LICENSE). Wflow is also available as +a [compiled executable](https://download.deltares.nl/en/download/wflow/) under the Deltares +terms and conditions. The wflow computational engine is built in the +[Julia](https://julialang.org/) language, a high-performance computing language. Wflow does not +include a graphical user interface and is designed for maximum user flexibility. Prior to 2021, +wflow was developed in Python on top of the PCRaster Python extension. The Python version is +[still available](https://github.com/openstreams/wflow), but not actively developed. ## Quick overview diff --git a/docs/model_docs/index.qmd b/docs/model_docs/index.qmd index c66ac2cf1..54ddc17ab 100644 --- a/docs/model_docs/index.qmd +++ b/docs/model_docs/index.qmd @@ -3,16 +3,16 @@ title: About --- As opposed to the user guide, which describes the steps needed to build and apply a model in -the software, this section explains the different model concepts that are available within -the modelling framework of wflow. Descriptions are given regarding the model concepts with -links to the original scientific papers which explain the concepts in more detail. The model -parameters which influence the processes are also shown, using inline code blocks. An -overview of all model parameters is also provided for easy reference, including their short -names, long descriptions and their units (see [parameters vertical concepts](@ref -params_vert) and [parameters lateral concepts](@ref params_lat)). +the software, this section explains the different model concepts that are available within the +modelling framework of wflow. Descriptions are given regarding the model concepts with links to +the original scientific papers which explain the concepts in more detail. The model parameters +which influence the processes are also shown, using inline code blocks. An overview of all +model parameters is also provided for easy reference, including their short names, long +descriptions and their units (see [parameters vertical concepts](@ref params_vert) and +[parameters lateral concepts](@ref params_lat)). ## Division between vertical and lateral -In the documentation we talk of `vertical` and `lateral` concepts. These are components in -the model that describe the vertical movement of water in each model grid cell and the -lateral movement of water across grid cells. +In the documentation we talk of `vertical` and `lateral` concepts. These are components in the +model that describe the vertical movement of water in each model grid cell and the lateral +movement of water across grid cells. diff --git a/docs/model_docs/lateral/gwf.qmd b/docs/model_docs/lateral/gwf.qmd index 33f4a8b6a..7e17e8f41 100644 --- a/docs/model_docs/lateral/gwf.qmd +++ b/docs/model_docs/lateral/gwf.qmd @@ -2,7 +2,8 @@ title: Groundwater flow --- -Single layer groundwater flow requires the four following components, and each is described in more detail below: +Single layer groundwater flow requires the four following components, and each is described in +more detail below: + aquifer + connectivity @@ -11,23 +12,23 @@ Single layer groundwater flow requires the four following components, and each i ## Aquifer types Groundwater flow can occur either in a confined or unconfined aquifer. Confined aquifers are -overlain by a poorly permeable confining layer (e.g. clay). No air can get in to fill the -pore space so that the aquifer always remains fully saturated. For a confined aquifer, water -will always flow along the complete height $\SIb{H}{m}$ over the aquifer and transmissivity +overlain by a poorly permeable confining layer (e.g. clay). No air can get in to fill the pore +space so that the aquifer always remains fully saturated. For a confined aquifer, water will +always flow along the complete height $\SIb{H}{m}$ over the aquifer and transmissivity $\SIb{kH}{m^2 d^{-1}}$ is a constant ($\SIb{k}{m d^{-1}}$ is the horizontal hydraulic conductivity). Specific storage is the amount of water an aquifer releases per unit change in hydraulic head, per unit volume of aquifer, as the aquifer and the groundwater itself is -compressed. Its value is much smaller than specific yield, between $10^{-5}$ (stiff) and $10^{-2}$ -(weak). +compressed. Its value is much smaller than specific yield, between $10^{-5}$ (stiff) and +$10^{-2}$ (weak). -The upper boundary of an unconfined aquifer is the water table (the phreatic -surface). Specific yield (or drainable porosity) represents the volumetric fraction the -aquifer will yield when all water drains and the pore volume is filled by air instead. -Specific yield will vary roughly between $0.05$ (clay) and $0.45$ (peat) (Johnson, 1967). +The upper boundary of an unconfined aquifer is the water table (the phreatic surface). Specific +yield (or drainable porosity) represents the volumetric fraction the aquifer will yield when +all water drains and the pore volume is filled by air instead. Specific yield will vary roughly +between $0.05$ (clay) and $0.45$ (peat) (Johnson, 1967). Groundwater flow is solved forward in time and central in space. The vertically averaged -governing equation for an inhomogeneous and isotropic aquifer in one dimension can be -written as: +governing equation for an inhomogeneous and isotropic aquifer in one dimension can be written +as: $$ S \frac{\partial \phi}{\partial t} = \frac{\partial}{\partial x} \left(kH \frac{\phi}{\delta x}\right) + Q @@ -35,8 +36,9 @@ $$ where $\SIb{S}{m m^{-1}}$ is storativity (or specific yield), $\SIb{\phi}{m}$ is hydraulic head, $t$ is time, $\SIb{k}{m t^{-1}}$ is horizontal hydraulic conductivity, $\SIb{H}{m}$ is -the (saturated) aquifer height: groundwater level - aquifer bottom elevation and $\SIb{Q}{m t^{-1}}$ represents fluxes from boundary conditions (e.g. recharge or abstraction), see -also [Aquifer boundary conditions](@ref). +the (saturated) aquifer height: groundwater level - aquifer bottom elevation and $\SIb{Q}{m +t^{-1}}$ represents fluxes from boundary conditions (e.g. recharge or abstraction), see also +[Aquifer boundary conditions](@ref). The simplest finite difference formulation is forward in time, central in space, and can be written as: @@ -45,10 +47,10 @@ $$ S_i \frac{\phi_{i}^{t+1} - \phi_i^{t}}{\Delta t} = -C_{i-1} (\phi_{i-1} - \phi_i) - C_i (\phi_{i+1} - \phi_i) + Q_i $$ -where $i$ is the cell index, $t$ is time, $\Delta t$ is the step size, $C_{i-1}$ -is the the intercell conductance between cell $i-1$ and $i$ and $C_i$ is the intercell -conductance between cell $i$ and $i+1$. The connection data between cells is stored as -part of the `Connectivity` struct, see also [Connectivity](@ref) for more information. +where $i$ is the cell index, $t$ is time, $\Delta t$ is the step size, $C_{i-1}$ is the the +intercell conductance between cell $i-1$ and $i$ and $C_i$ is the intercell conductance between +cell $i$ and $i+1$. The connection data between cells is stored as part of the `Connectivity` +struct, see also [Connectivity](@ref) for more information. Conductance $C$ is defined as: @@ -56,30 +58,29 @@ $$ C = \frac{kH w}{l} $$ -where $\SIb{w}{m}$ is the width of the cell to cell connection, and $\SIb{l}{m}$ is the length of -the cell to cell connection. $k$ and $H$ may both vary in space; intercell conductance -is therefore an average using the properties of two cells. For the calculation of the -intercell conductance $C$ the harmonic mean is used (see also Goode and Appel, 1992), here -between cell index $i$ and cell index $i+1$, in the $x$ direction: +where $\SIb{w}{m}$ is the width of the cell to cell connection, and $\SIb{l}{m}$ is the length +of the cell to cell connection. $k$ and $H$ may both vary in space; intercell conductance is +therefore an average using the properties of two cells. For the calculation of the intercell +conductance $C$ the harmonic mean is used (see also Goode and Appel, 1992), here between cell +index $i$ and cell index $i+1$, in the $x$ direction: $$ C_i = w \frac{k_iH_i\cdot k_{i+1}H_{i+1}}{k_iH_i \cdot l_{i+1} + k_{i+1}H_{i+1} \cdot l_i} $$ -where $\SIb{H}{m}$ is the aquifer top - aquifer bottom, and $k$, $l_i$ is the length in -cell $i$ ($0.5 \Delta x_i$), $l_{i+1}$ is the length in cell $i+1$ ($0.5 \Delta -x_{i+1}$) and $w$ as previously defined. For an unconfined aquifer the intercell -conductance is scaled by using the "upstream saturated fraction" as the MODFLOW -documentation calls it. In this approach, the saturated thickness of a cell-to-cell is -approximated using the cell with the highest head. This results in a consistent -overestimation of the saturated thickness, but it avoids complexities related with cell -drying and rewetting, such as having to define a "wetting threshold" or a "wetting factor". -See also the documentation for MODFLOW-NWT (Niswonger et al., 2011) or MODFLOW6 (Langevin et -al., 2017) for more background information. For more background on drying and rewetting, see -for example McDonald et al. (1991). - -For the finite difference formulation, there is only one unknown, $\phi_i^{t+1}$. -Reshuffling terms: +where $\SIb{H}{m}$ is the aquifer top - aquifer bottom, and $k$, $l_i$ is the length in cell +$i$ ($0.5 \Delta x_i$), $l_{i+1}$ is the length in cell $i+1$ ($0.5 \Delta x_{i+1}$) and $w$ +as previously defined. For an unconfined aquifer the intercell conductance is scaled by using +the "upstream saturated fraction" as the MODFLOW documentation calls it. In this approach, the +saturated thickness of a cell-to-cell is approximated using the cell with the highest head. +This results in a consistent overestimation of the saturated thickness, but it avoids +complexities related with cell drying and rewetting, such as having to define a "wetting +threshold" or a "wetting factor". See also the documentation for MODFLOW-NWT (Niswonger et al., +2011) or MODFLOW6 (Langevin et al., 2017) for more background information. For more background +on drying and rewetting, see for example McDonald et al. (1991). + +For the finite difference formulation, there is only one unknown, $\phi_i^{t+1}$. Reshuffling +terms: $$ \phi_i^{t+1} = \phi_i^t + (C_{i-1} (\phi_i - \phi_{i-1}) + C_i (\phi_{i+1} - \phi_i) + Q_i) \frac{Δt}{S_i} @@ -93,12 +94,12 @@ $$ \frac{\Delta t k H}{\Delta x \Delta y S} \le \frac{1}{4} $$ -where $\SIb{\Delta t}{d}$ is the stable time step size, $\SIb{\Delta x}{m}$ is the cell length in -the $x$ direction and $\SIb{\Delta y}{m}$ is the cell length in the $y$ direction, $\SIb{k}{m^2 d^{-1}}$ is -the horizontal hydraulic conductivity and $\SIb{H}{m}$ is the saturated -thickness of the aquifer. For each cell $\frac{\Delta x \Delta y S}{k H}$ is -calculated, the minimum of these values is determined, and multiplied by $\frac{1}{4}$, to -get the stable time step size. +where $\SIb{\Delta t}{d}$ is the stable time step size, $\SIb{\Delta x}{m}$ is the cell length +in the $x$ direction and $\SIb{\Delta y}{m}$ is the cell length in the $y$ direction, +$\SIb{k}{m^2 d^{-1}}$ is the horizontal hydraulic conductivity and $\SIb{H}{m}$ is the +saturated thickness of the aquifer. For each cell $\frac{\Delta x \Delta y S}{k H}$ is +calculated, the minimum of these values is determined, and multiplied by $\frac{1}{4}$, to get +the stable time step size. For more details about the finite difference formulation and the stable time step size criterion we refer to the paper of Chu and Willis (1984). @@ -132,8 +133,8 @@ Dirichlet boundary conditions can be specified through the field `constanthead` end ``` -For the model `SBM + Groundwater flow` this boundary condition is optional, and if used -should be specified in the TOML file as follows (see also +For the model `SBM + Groundwater flow` this boundary condition is optional, and if used should +be specified in the TOML file as follows (see also [sbm\_gwf\_config.toml](https://github.com/Deltares/Wflow.jl/blob/master/test/sbm_gwf_config.toml)): ```toml @@ -148,7 +149,7 @@ The flux between river and aquifer is calculated using Darcy's law following the MODFLOW: $$ - \subtext{Q}{riv} = + \subtext{Q}{riv} = \begin{align*} \begin{cases} C_i \min \left\{\subtext{h}{riv} - \subtext{B}{riv}, \subtext{h}{riv} - \phi\right\} &\text{ if }\quad \subtext{h}{riv} > \phi \\ @@ -157,16 +158,18 @@ $$ \end{align*} $$ -where $\SIb{\subtext{Q}{riv}}{L^3 T^{-1}}$ is the exchange flux from river to aquifer, $\SIb{C_i}{L^2 T^{-1}}$ is the river bed infiltration conductance, $\SIb{C_e}{L^2 T^{-1}}$ -is the river bed exfiltration conductance, $\SIb{\subtext{B}{riv}}{L}$ the bottom of the river bed, -$\SIb{\subtext{h}{riv}}{L}$ is the river stage and $\SIb{\phi}{L}$ is the hydraulic head in the river cell. +where $\SIb{\subtext{Q}{riv}}{L^3 T^{-1}}$ is the exchange flux from river to aquifer, +$\SIb{C_i}{L^2 T^{-1}}$ is the river bed infiltration conductance, $\SIb{C_e}{L^2 T^{-1}}$ is +the river bed exfiltration conductance, $\SIb{\subtext{B}{riv}}{L}$ the bottom of the river +bed, $\SIb{\subtext{h}{riv}}{L}$ is the river stage and $\SIb{\phi}{L}$ is the hydraulic head +in the river cell. -The Table in the Groundwater flow [river boundary condition](@ref gwf_river_params) section -of the Model parameters provides the parameters of the struct `River`. Parameters that can -be set directly from the static input data (netCDF) are marked in this Table. +The Table in the Groundwater flow [river boundary condition](@ref gwf_river_params) section of +the Model parameters provides the parameters of the struct `River`. Parameters that can be set +directly from the static input data (netCDF) are marked in this Table. -The exchange flux (river to aquifer) $\subtext{Q}{riv}$ is an output variable (field `flux` of the -`River` struct), and is used to update the total flux in a river cell. For the model `SBM + +The exchange flux (river to aquifer) $\subtext{Q}{riv}$ is an output variable (field `flux` of +the `River` struct), and is used to update the total flux in a river cell. For the model `SBM + Groundwater flow`, the water level $\SIb{h}{m}$ of the river kinematic wave in combination with the river `bottom` is used to update the `stage` field of the `River` struct each time step. @@ -178,15 +181,15 @@ $$ $$ where $\SIb{\subtext{Q}{drain}}{L^3 T^{-1}}$ is the exchange flux from drains to aquifer, -$\SIb{\subtext{C}{drain}}{L^2 T^{-1}}$ is the drain conductance, $\SIb{\subtext{h}{drain}}{L}$ is the drain -elevation and $\SIb{\phi}{L}$ is the hydraulic head in the cell with drainage. +$\SIb{\subtext{C}{drain}}{L^2 T^{-1}}$ is the drain conductance, $\SIb{\subtext{h}{drain}}{L}$ +is the drain elevation and $\SIb{\phi}{L}$ is the hydraulic head in the cell with drainage. The table in the Groundwater flow [drainage boundary condition](@ref gwf_drainage_params) section of the Model parameters provides the parameters of the struct `Drainage`. Parameters that can be set directly from the static input data (netCDF) are marked in this Table. -The exchange flux (drains to aquifer) $\subtext{Q}{drain}$ is an output variable (field `flux` of -struct `Drainage`), and is used to update the total flux in a cell with drains. For the +The exchange flux (drains to aquifer) $\subtext{Q}{drain}$ is an output variable (field `flux` +of struct `Drainage`), and is used to update the total flux in a cell with drains. For the model `SBM + Groundwater flow` this boundary condition is optional, and if used should be specified in the TOML file as follows (see also [sbm\_gwf\_config.toml](https://github.com/Deltares/Wflow.jl/blob/master/test/sbm_gwf_config.toml)): @@ -203,25 +206,24 @@ $$ Q_r = R \, A $$ -with $\SIb{}{L T^{-1}}$ the recharge rate and $\SIb{A}{L^2}$ the area of the aquifer -cell. +with $\SIb{}{L T^{-1}}$ the recharge rate and $\SIb{A}{L^2}$ the area of the aquifer cell. The table in the Groundwater flow [recharge boundary condition](@ref gwf_recharge_params) section of the Model parameters section provides the parameters of the struct `Recharge`. Parameters that can be set directly from the static input data (netCDF) are marked in this Table. -The recharge flux $Q_r$ is an output variable (field `flux` of struct `Recharge`), and is -used to update the total flux in a cell where recharge occurs. For the model `SBM + -Groundwater flow`, the recharge rate from the vertical SBM concept `recharge` [mm] is used -to update the `rate` field of the `Recharge` struct each time step. The `rate` field is -multiplied by the `area` field of the aquifer. +The recharge flux $Q_r$ is an output variable (field `flux` of struct `Recharge`), and is used +to update the total flux in a cell where recharge occurs. For the model `SBM + Groundwater +flow`, the recharge rate from the vertical SBM concept `recharge` [mm] is used to update the +`rate` field of the `Recharge` struct each time step. The `rate` field is multiplied by the +`area` field of the aquifer. ### Head boundary -This boundary is a fixed head with time (not affected by the model stresses over time)) -outside of the model domain, and is generally used to avoid an unnecessary extension of the -model domain to the location of the fixed boundary (for example a large lake). The flux from -the boundary $\SIb{Q_{hb}}{L^3 T^{-1}}$ is calculated as follows: +This boundary is a fixed head with time (not affected by the model stresses over time) outside +of the model domain, and is generally used to avoid an unnecessary extension of the model +domain to the location of the fixed boundary (for example a large lake). The flux from the +boundary $\SIb{Q_{hb}}{L^3 T^{-1}}$ is calculated as follows: $$ Q_{hb} = C_{hb} (\phi_{hb} - \phi) @@ -233,10 +235,9 @@ head of the head boundary and $\phi$ the head of the aquifer cell. The table in the Groundwater flow [head boundary condition](@ref gwf_headboundary_params) section of the Model parameters provides the parameters of the struct `HeadBoundary`. -The head boundary flux $Q_{hb}$ is an output variable (field `flux` of struct -`HeadBoundary`), and is used to update the total flux in a cell where this type of boundary -occurs. The parameter Head $\phi_{hb}$ can be specified as a fixed or time dependent -value. +The head boundary flux $Q_{hb}$ is an output variable (field `flux` of struct `HeadBoundary`), +and is used to update the total flux in a cell where this type of boundary occurs. The +parameter Head $\phi_{hb}$ can be specified as a fixed or time dependent value. ::: {.callout-note} This boundary is not (yet) part of the model `SBM + Groundwater flow`. @@ -248,34 +249,33 @@ A volumetric well rate $\SIb{}{L^3 T^{-1}}$ can be specified as a boundary condi The Table in the [well boundary condition](@ref well_boundary_params) section of the Model parameters provides the parameters of the struct `Well`. -The volumetric well rate $\subtext{Q}{well}$ can be can be specified as a fixed or time dependent -value. If a cell is dry, the actual well flux `flux` is set to zero (see also the last note -on this page). +The volumetric well rate $\subtext{Q}{well}$ can be can be specified as a fixed or time +dependent value. If a cell is dry, the actual well flux `flux` is set to zero (see also the +last note on this page). ::: {.callout-note} This boundary is not (yet) part of the model `SBM + Groundwater flow`. ::: ::: {.callout-note} -For an unconfined aquifer the boundary fluxes are checked, in case of a dry aquifer cell -a negative flux is not allowed. +For an unconfined aquifer the boundary fluxes are checked, in case of a dry aquifer cell a +negative flux is not allowed. ::: ## References -+ Chu, W. S., & Willis, R. (1984). An explicit finite difference model for unconfined - aquifers. Groundwater, 22(6), 728-734. ++ Chu, W. S., & Willis, R. (1984). An explicit finite difference model for unconfined aquifers. + Groundwater, 22(6), 728-734. + Goode, D. J., & Appel, C. A. (1992). Finite-Difference Interblock Transmissivity for Unconfined Aquifers and for Aquifers having Smoothly Varying Transmissivity Water-resources investigations report, 92, 4124. -+ Johnson, A. I. (1967), Specific yield: compilation of specific yields for various - materials, Water Supply Paper 1662-D, Washington, D.C.: U.S. Government Printing Office, - p. 74, doi:10.3133/wsp1662D. -+ Langevin, C.D., Hughes, J.D., Banta, E.R., Niswonger, R.G., Panday, Sorab, and Provost, - A.M., 2017, Documentation for the MODFLOW 6 Groundwater Flow Model: U.S. Geological Survey ++ Johnson, A. I. (1967), Specific yield: compilation of specific yields for various materials, + Water Supply Paper 1662-D, Washington, D.C.: U.S. Government Printing Office, p. 74, + doi:10.3133/wsp1662D. ++ Langevin, C.D., Hughes, J.D., Banta, E.R., Niswonger, R.G., Panday, Sorab, and Provost, A.M., + 2017, Documentation for the MODFLOW 6 Groundwater Flow Model: U.S. Geological Survey Techniques and Methods, book 6, chap. A55, 197 p., https://doi.org/10.3133/tm6A55. -+ McDonald, M.G., Harbaugh, A.W., Orr, B.R., and Ackerman, D.J., 1991, A method of - converting no-flow cells to variable-head cells for the U.S. Geological Survey modular - finite-difference groundwater flow model: U.S. Geological Survey Open-File Report 91-536, - 99 p. -+ Niswonger, R.G., Panday, Sorab, and Ibaraki, Motomu, 2011, MODFLOW-NWT, A Newton - formulation for MODFLOW-2005: U.S. Geological Survey Techniques and Methods 6-A37, 44 p. \ No newline at end of file ++ McDonald, M.G., Harbaugh, A.W., Orr, B.R., and Ackerman, D.J., 1991, A method of converting + no-flow cells to variable-head cells for the U.S. Geological Survey modular finite-difference + groundwater flow model: U.S. Geological Survey Open-File Report 91-536, 99 p. ++ Niswonger, R.G., Panday, Sorab, and Ibaraki, Motomu, 2011, MODFLOW-NWT, A Newton formulation + for MODFLOW-2005: U.S. Geological Survey Techniques and Methods 6-A37, 44 p. \ No newline at end of file diff --git a/docs/model_docs/lateral/kinwave.qmd b/docs/model_docs/lateral/kinwave.qmd index defd201dd..0c3e4ccb5 100644 --- a/docs/model_docs/lateral/kinwave.qmd +++ b/docs/model_docs/lateral/kinwave.qmd @@ -3,9 +3,9 @@ title: Kinematic wave --- ## Surface routing -The main flow routing scheme available in `Wflow.jl` is the kinematic wave approach for -channel and overland flow, assuming that water flow is mostly controlled by topography. The -kinematic wave equations are (Chow, 1988): +The main flow routing scheme available in `Wflow.jl` is the kinematic wave approach for channel +and overland flow, assuming that water flow is mostly controlled by topography. The kinematic +wave equations are (Chow, 1988): $$ \begin{gathered} @@ -20,17 +20,16 @@ $$ \dfrac{\partial Q}{\partial x} + \alpha \beta Q^{\beta - 1} \dfrac{\partial Q}{\partial t} = q. $$ -Here $\SIb{Q}{m^3 s^{-1}}$ is the surface runoff in the kinematic wave, $\SIb{x}{m}$ is the length of -the runoff pathway, $\SIb{A}{m}$ is the cross-section area of the runoff pathway, +Here $\SIb{Q}{m^3 s^{-1}}$ is the surface runoff in the kinematic wave, $\SIb{x}{m}$ is the +length of the runoff pathway, $\SIb{A}{m}$ is the cross-section area of the runoff pathway, $\SIb{t}{s}$ is the integration timestep and $\alpha$ and $\beta$ are unitless coefficients. These equations are solved with a nonlinear scheme using Newton's method and can also be iterated depending on the model space and time resolution. By default, the iterations are -performed until a stable solution is reached ($\epsilon < 10^{-12}$). For larger models, -the number of iterations can also be fixed for to a specific sub-timestep (in seconds) for -both overland and channel flows to improve simulation time. To enable (fixed or not) -iterations of the kinematic wave the following lines can be inserted in the TOML file of the -model: +performed until a stable solution is reached ($\epsilon < 10^{-12}$). For larger models, the +number of iterations can also be fixed for to a specific sub-timestep (in seconds) for both +overland and channel flows to improve simulation time. To enable (fixed or not) iterations of +the kinematic wave the following lines can be inserted in the TOML file of the model: ```toml [model] @@ -42,9 +41,9 @@ kw_river_tstep = 900 kw_land_tstep = 3600 ``` -The $\alpha$ parameter of the kinematic wave is fixed. To estimate the wetted perimeter -for the calculation of the $\alpha$ parameter a bankfull river depth map (default value -is $\SI{1.0}{m}$) for the river can be provided as follows: +The $\alpha$ parameter of the kinematic wave is fixed. To estimate the wetted perimeter for the +calculation of the $\alpha$ parameter a bankfull river depth map (default value is +$\SI{1.0}{m}$) for the river can be provided as follows: ```toml [input.lateral.river] @@ -59,54 +58,54 @@ river kinematic wave network. ## Inflow External water (supply/abstraction) `inflow` $\SIb{}{m^3 s^{-1}}$ can be added to the -kinematic wave for surface water routing, as a cyclic parameter or as part of forcing (see -also [Input section](@ref)). +kinematic wave for surface water routing, as a cyclic parameter or as part of forcing (see also +[Input section](@ref)). ## Abstractions Abstractions from the river through the variable `abstraction` $\SIb{}{m^3 s^{-1}}$ are -possible when water demand and allocation is computed. The variable `abstraction` is set -from the water demand and allocation module each time step. The `abstraction` is divided by -the length of the runoff pathway and subtracted from the lateral inflow of the kinematic -wave routing scheme for river flow. +possible when water demand and allocation is computed. The variable `abstraction` is set from +the water demand and allocation module each time step. The `abstraction` is divided by the +length of the runoff pathway and subtracted from the lateral inflow of the kinematic wave +routing scheme for river flow. ## Subsurface flow routing In the SBM model the kinematic wave approach is used to route subsurface flow laterally. -Different vertical hydraulic conductivity depth profiles are possible as part of the -vertical [SBM](@ref soil) concept, and these profiles (after unit conversion) are also used -to compute lateral subsurface flow. The following profiles (see [SBM](@ref soil) for a -detailed description) are available: +Different vertical hydraulic conductivity depth profiles are possible as part of the vertical +[SBM](@ref soil) concept, and these profiles (after unit conversion) are also used to compute +lateral subsurface flow. The following profiles (see [SBM](@ref soil) for a detailed +description) are available: - `exponential` (default) - `exponential_constant` - `layered` -- `layered_exponential` +- `layered_exponential` -For the profiles `exponential` and `exponential_constant`, the saturated store ``S`` is -drained laterally by saturated downslope subsurface flow for a slope with width ``\SIb{w}{m}`` +For the profiles `exponential` and `exponential_constant`, the saturated store ``S`` is drained +laterally by saturated downslope subsurface flow for a slope with width ``\SIb{w}{m}`` according to: $$ Q = K_0\tan(\beta)w\begin{cases} - \frac{1}{f}\left(e^{-fz_i}-e^{-f\subtext{z}{exp}}\right) + + \frac{1}{f}\left(e^{-fz_i}-e^{-f\subtext{z}{exp}}\right) + e^{-f\subtext{z}{exp}}(z_t-\subtext{z}{exp}) & \text{if $z_i < \subtext{z}{exp}$}\\ \\ e^{-f\subtext{z}{exp}}(z_t - z_i) & \text{if $z_i \ge \subtext{z}{exp}$}, \end{cases} $$ -where $\beta$ is element slope angle, $\SIb{Q}{m^3 d^{-1}}$ is subsurface flow, -$\SIb{K_0}{m d^{-1}}$ is the saturated hydraulic conductivity at the soil surface, $\SIb{z_i}{m}$ -is the water table depth, $\SIb{z_{t}}{m}$ is the total soil depth, $\SIb{f}{m^{-1}}$ is a scaling -parameter that controls the decrease of $K_0$ with depth and -$\SIb{\subtext{z}{exp}}{m}$ is the depth from soil surface for which the exponential decline of -$K_0$ is valid. For the `exponential` profile, $\subtext{z}{exp}$ is equal to $z_t$. +where $\beta$ is element slope angle, $\SIb{Q}{m^3 d^{-1}}$ is subsurface flow, $\SIb{K_0}{m +d^{-1}}$ is the saturated hydraulic conductivity at the soil surface, $\SIb{z_i}{m}$ is the +water table depth, $\SIb{z_{t}}{m}$ is the total soil depth, $\SIb{f}{m^{-1}}$ is a scaling +parameter that controls the decrease of $K_0$ with depth and $\SIb{\subtext{z}{exp}}{m}$ is the +depth from soil surface for which the exponential decline of $K_0$ is valid. For the +`exponential` profile, $\subtext{z}{exp}$ is equal to $z_t$. Combining with the following continuity equation: $$ (\theta_s-\theta_r)w\frac{\partial h}{\partial t} = -\frac{\partial Q}{\partial x} + wr $$ -where $\SIb{h}{m}$ is the water table height, $\SIb{x}{m}$ is the distance downslope, and $\SIb{r}{m d^{-1}}$ is -the net input rate to the saturated store. Substituting for $h +where $\SIb{h}{m}$ is the water table height, $\SIb{x}{m}$ is the distance downslope, and +$\SIb{r}{m d^{-1}}$ is the net input rate to the saturated store. Substituting for $h (\frac{\partial Q}{\partial h})$, gives: $$ @@ -124,7 +123,8 @@ $$ $$ For the `layered` and `layered_exponential` profiles the equivalent horizontal hydraulic -conductivity $\SIb{K_h}{m d^{-1}}$ is calculated for water table height $\SIb{h = z_t-z_i}{m}$ and lateral subsurface flow is calculated as follows: +conductivity $\SIb{K_h}{m d^{-1}}$ is calculated for water table height $\SIb{h = z_t-z_i}{m}$ +and lateral subsurface flow is calculated as follows: $$ Q = K_h h \tan(\beta) w, @@ -140,18 +140,17 @@ The kinematic wave equation for lateral subsurface flow is solved iteratively us method. ::: {.callout-note} - For the lateral subsurface flow kinematic wave the model timestep is not adjusted. - For certain model timestep and model grid size combinations this may result in loss of - accuracy. +For the lateral subsurface flow kinematic wave the model timestep is not adjusted. For certain +model timestep and model grid size combinations this may result in loss of accuracy. ::: ## Multi-Threading The kinematic wave calculations for surface - and subsurface flow routing can be executed in -parallel using multiple threads. In the model section of the TOML file, a minimum stream -order can be provided to define subbasins for the river (default is 6) and land domain -(default is 5). Subbasins are created at all confluences where each branch has a minimal -stream order. Based on the subbasins a directed acyclic graph is created that controls the -order of execution and which subbasins can run in parallel. +parallel using multiple threads. In the model section of the TOML file, a minimum stream order +can be provided to define subbasins for the river (default is 6) and land domain (default is +5). Subbasins are created at all confluences where each branch has a minimal stream order. +Based on the subbasins a directed acyclic graph is created that controls the order of execution +and which subbasins can run in parallel. ```toml [model] @@ -160,11 +159,11 @@ min_streamorder_land = 4 # minimum stream order to delineate subbasins for land ``` ## Subcatchment flow -Normally the the kinematic wave is continuous throughout the model. By using the `pits` -entry in the model and input sections of the TOML file all flow is at the subcatchment only -(upstream of the pit locations, defined by the netCDF variable `wflow_pits` in the example -below) and no flow is transferred from one subcatchment to another. This can be convenient -when connecting the result of the model to a water allocation model such as Ribasim. +Normally the the kinematic wave is continuous throughout the model. By using the `pits` entry +in the model and input sections of the TOML file all flow is at the subcatchment only (upstream +of the pit locations, defined by the netCDF variable `wflow_pits` in the example below) and no +flow is transferred from one subcatchment to another. This can be convenient when connecting +the result of the model to a water allocation model such as Ribasim. ```toml [input] @@ -176,25 +175,25 @@ pits = true ``` ## Limitations -The kinematic wave approach for channel, overland and lateral subsurface flow, assumes that -the topography controls water flow mostly. This assumption holds for steep terrain, but in -less steep terrain the hydraulic gradient is likely not equal to the surface slope -(subsurface flow), or pressure differences and inertial momentum cannot be neglected -(channel and overland flow). In addition, while the kinematic wave equations are solved -with a nonlinear scheme using Newton's method (Chow, 1988), other model equations are solved -through a simple explicit scheme. In summary the following limitations apply: - -+ Channel flow, and to a lesser degree overland flow, may be unrealistic in terrain that is - not steep, and where pressure forces and inertial momentum are important. +The kinematic wave approach for channel, overland and lateral subsurface flow, assumes that the +topography controls water flow mostly. This assumption holds for steep terrain, but in less +steep terrain the hydraulic gradient is likely not equal to the surface slope (subsurface +flow), or pressure differences and inertial momentum cannot be neglected (channel and overland +flow). In addition, while the kinematic wave equations are solved with a nonlinear scheme using +Newton's method (Chow, 1988), other model equations are solved through a simple explicit +scheme. In summary the following limitations apply: + ++ Channel flow, and to a lesser degree overland flow, may be unrealistic in terrain that is not + steep, and where pressure forces and inertial momentum are important. + The lateral movement of subsurface flow may be very wrong in terrain that is not steep. ## External inflows -External inflows, for example water supply or abstractions, can be added to the kinematic -wave via the `inflow` variable. For this, the user can supply a 2D map of the inflow, as a -cyclic parameter or as part of forcing (see also [Input section](@ref)). These inflows are -added or abstracted from the upstream inflow `qin` before running the kinematic wave to -solve the impact on resulting `q`. In case of a negative inflow (abstractions), a minimum of -zero is applied to the upstream flow `qin`. +External inflows, for example water supply or abstractions, can be added to the kinematic wave +via the `inflow` variable. For this, the user can supply a 2D map of the inflow, as a cyclic +parameter or as part of forcing (see also [Input section](@ref)). These inflows are added or +abstracted from the upstream inflow `qin` before running the kinematic wave to solve the impact +on resulting `q`. In case of a negative inflow (abstractions), a minimum of zero is applied to +the upstream flow `qin`. ## References + Chow, V., Maidment, D. and Mays, L., 1988, Applied Hydrology. McGraw-Hill Book Company, diff --git a/docs/model_docs/lateral/local-inertial.qmd b/docs/model_docs/lateral/local-inertial.qmd index da9fa0c82..8f9bc6466 100644 --- a/docs/model_docs/lateral/local-inertial.qmd +++ b/docs/model_docs/lateral/local-inertial.qmd @@ -4,18 +4,19 @@ title: Local inertial ## River and floodplain routing The local inertial approximation of shallow water flow neglects only the convective -acceleration term in the Saint-Venant momentum conservation equation. The numerical solution -of the local inertial approximation on a staggered grid is as follows (Bates et al., 2010): +acceleration term in the Saint-Venant momentum conservation equation. The numerical solution of +the local inertial approximation on a staggered grid is as follows (Bates et al., 2010): $$ Q_{t+\Delta t} = \frac{Q_t - g A_t \Delta t S_t}{(1+g\Delta t n^2 |Q_t| / (R_t^{4/3} A_t))} $$ -where $\SIb{Q_{t+\Delta t}}{m^3 s^{-1}}$ is the river flow at time step $t+\Delta t$, $\SIb{g}{m s^{-2}}$ is -acceleration due to gravity, $\SIb{A_t}{m^2}$ is the cross sectional flow area at the -previous time step, $\SIb{R_t}{m}$ is the hydraulic radius at the previous time step, $\SIb{Q_t}{m^3 s^{-1}}$ is -the river flow at the previous time step, $S_t$ is the water surface slope at -the previous time step and $\SIb{n}{m^{-\frac{1}{3}} s}$ is the Manning's roughness coefficient. +where $\SIb{Q_{t+\Delta t}}{m^3 s^{-1}}$ is the river flow at time step $t+\Delta t$, +$\SIb{g}{m s^{-2}}$ is acceleration due to gravity, $\SIb{A_t}{m^2}$ is the cross sectional +flow area at the previous time step, $\SIb{R_t}{m}$ is the hydraulic radius at the previous +time step, $\SIb{Q_t}{m^3 s^{-1}}$ is the river flow at the previous time step, $S_t$ is the +water surface slope at the previous time step and $\SIb{n}{m^{-\frac{1}{3}} s}$ is the +Manning's roughness coefficient. The momentum equation is applied to each link between two river grid cells, while the continuity equation over $\Delta t$ is applied to each river cell: @@ -24,10 +25,10 @@ $$ h^{t+\Delta t} = h^t + \Delta t \frac{\subtext{Q^{t+\Delta t}}{src} - \subtext{Q^{t+\Delta t}}{dst}}{A} $$ -where $\SIb{h^{t+\Delta t}}{m}$ is the water depthat time step $t+\Delta t$, $\SIb{h^t}{m}$ is the -water depth at the previous time step, $\SIb{A}{m^2}$ is the river area and $\SIb{\subtext{Q}{src}}{m^3 s^{-1}}$ -and $\SIb{\subtext{Q}{dst}}{m^3 s^{-1}}$ represent river flow at the upstream and downstream link of the -river cell, respectively. +where $\SIb{h^{t+\Delta t}}{m}$ is the water depthat time step $t+\Delta t$, $\SIb{h^t}{m}$ is +the water depth at the previous time step, $\SIb{A}{m^2}$ is the river area and +$\SIb{\subtext{Q}{src}}{m^3 s^{-1}}$ and $\SIb{\subtext{Q}{dst}}{m^3 s^{-1}}$ represent river +flow at the upstream and downstream link of the river cell, respectively. The model time step $\Delta t$ for the local inertial model is estimated based on the Courant-Friedrichs-Lewy condition (Bates et al., 2010): @@ -37,11 +38,11 @@ $$ $$ where $\sqrt{gh_i}$ is the wave celerity for river cell $i$ , $\SIb{\Delta x_i}{m}$ is the -river length for river cell $i$ and $\alpha$ is a coefficient (typically between $0.2$ -and $0.7$) to enhance the stability of the simulation. +river length for river cell $i$ and $\alpha$ is a coefficient (typically between $0.2$ and +$0.7$) to enhance the stability of the simulation. -In the TOML file the following properties related to the local inertial model can be -provided for the `sbm` and `sbm_gwf` model types: +In the TOML file the following properties related to the local inertial model can be provided +for the `sbm` and `sbm_gwf` model types: ```toml [model] @@ -52,59 +53,59 @@ h_thresh = 0.1 # water depth [m] threshold for calculating fl floodplain_1d = true # include 1D floodplain schematization (default = false) ``` -Two optional constant boundary conditions `riverlength_bc` and `riverdepth_bc` can be -provided at a river outlet node (or multiple river outlet nodes) through the model parameter -netCDF file, as follows: +Two optional constant boundary conditions `riverlength_bc` and `riverdepth_bc` can be provided +at a river outlet node (or multiple river outlet nodes) through the model parameter netCDF +file, as follows: ```toml [input.lateral.river] riverlength_bc = "riverlength_bc" # optional river length [m], default = 1e04 riverdepth_bc = "riverdepth_bc" # optional river depth [m], default = 0.0 ``` -These boundary conditions are copied to a ghost node (downstream of the river outlet node) -in the code. +These boundary conditions are copied to a ghost node (downstream of the river outlet node) in +the code. -The optional 1D floodplain schematization is based on provided flood volumes as a function -of flood depth (per flood depth interval) for each river cell. Wflow calculates from these -flood volumes a rectangular floodplain profile for each flood depth interval. Routing is -done separately for the river channel and floodplain. +The optional 1D floodplain schematization is based on provided flood volumes as a function of +flood depth (per flood depth interval) for each river cell. Wflow calculates from these flood +volumes a rectangular floodplain profile for each flood depth interval. Routing is done +separately for the river channel and floodplain. The momentum equation is most stable for low slope environments, and to keep the simulation stable for (partly) steep environments the `froude_limit` option is set to true by default. This setting limits flow conditions to subcritical-critical conditions based on the Froude -number ($\le 1$), similar to Coulthard et al. (2013) in the CAESAR-LISFLOOD model and Adams -et al. (2017) in the Landlab v1.0 OverlandFlow component. The froude number $\mathrm{Fr}$ on a link -is calculated as follows: +number ($\le 1$), similar to Coulthard et al. (2013) in the CAESAR-LISFLOOD model and Adams et +al. (2017) in the Landlab v1.0 OverlandFlow component. The froude number $\mathrm{Fr}$ on a +link is calculated as follows: $$ \mathrm{Fr} = \frac{u}{\sqrt{gh_f}} $$ -where $\sqrt{gh_f}$ is the wave celerity on a link and $u$ is the water velocity on a -link. If the water velocity from the local inertial model is causing the Froude number to be -greater than $1.0$ , the water velocity (and flow) is reduced in order to maintain a Froude -number of $1.0$. +where $\sqrt{gh_f}$ is the wave celerity on a link and $u$ is the water velocity on a link. If +the water velocity from the local inertial model is causing the Froude number to be greater +than $1.0$ , the water velocity (and flow) is reduced in order to maintain a Froude number of +$1.0$. -The downstream boundary condition basically simulates a zero water depth boundary condition -at a set distance, as follows. For the downstream boundary condition (ghost point) the river +The downstream boundary condition basically simulates a zero water depth boundary condition at +a set distance, as follows. For the downstream boundary condition (ghost point) the river width, river bed elevation and Manning's roughness coefficient are copied from the upstream -river cell. The river length $\SIb{}{m}$ of the boundary cell can be set through the TOML file with -`riverlength_bc`, and has a default value of $\SI{10}{km}$. The water depth at the boundary cell is -fixed at $\SI{0.0}{m}$. +river cell. The river length $\SIb{}{m}$ of the boundary cell can be set through the TOML file +with `riverlength_bc`, and has a default value of $\SI{10}{km}$. The water depth at the +boundary cell is fixed at $\SI{0.0}{m}$. Simplified [reservoir and lake](@ref reservoir_lake) models can be included as part of the local inertial model for river flow (1D) and river and overland flow combined (see next -section). Reservoir and lake models are included as a boundary point with zero water depth -for both river and overland flow. For river flow the reservoir or lake model replaces the -local inertial model at the reservoir or lake location, and $Q$ is set by the outflow from -the reservoir or lake. Overland flow at a reservoir or lake location is not allowed to or -from the downstream river grid cell. +section). Reservoir and lake models are included as a boundary point with zero water depth for +both river and overland flow. For river flow the reservoir or lake model replaces the local +inertial model at the reservoir or lake location, and $Q$ is set by the outflow from the +reservoir or lake. Overland flow at a reservoir or lake location is not allowed to or from the +downstream river grid cell. ## Overland flow (2D) -For the simulation of 2D overland flow on a staggered grid the numerical scheme proposed by -de Almeida et al. (2012) is adopted. The explicit solution for the estimation of water -discharge between two cells in the x-direction is of the following form (following the -notation of Almeida et al. (2012)): +For the simulation of 2D overland flow on a staggered grid the numerical scheme proposed by de +Almeida et al. (2012) is adopted. The explicit solution for the estimation of water discharge +between two cells in the x-direction is of the following form (following the notation of +Almeida et al. (2012)): $$ \begin{split} @@ -114,22 +115,23 @@ Q_{(i+1/2)}^{n})\right]- g h_f \frac{\Delta t}{\Delta x} (\eta^n_i - \eta^n_{i-1 $$ where subscripts $i$ and $n$ refer to space and time indices, respectively. Subscript -$i-\frac{1}{2}$ is to the link between node $i$ and $i-1$, subscript $i+\frac{1}{2}$ is the link -between node $i$ and node $i+1$, and subscript $i-\frac{3}{2}$ is the link between node $i-1$ -and node $i-2$. $\SIb{Q}{m^3 s^{-1}}$ is the water discharge, $\SIb{\eta}{m}$ is the water -surface elevation, $\SIb{h_f}{m}$ is the water depth between cells, $\SIb{n}{m^{-\frac{1}{3}} s}$ is the Manning's -roughness coefficient, $\SIb{g}{m s^{-2}}$ is acceleration due to gravity, -$\SIb{\Delta t}{s}$ is the adaptive model time step, $\SIb{\Delta x}{m}$ is the distance between -two cells and $\SIb{\Delta y}{m}$ is the flow width. Below the staggered grid and variables of -the numerical solution in the x-direction, based on Almeida et al. (2012): +$i-\frac{1}{2}$ is to the link between node $i$ and $i-1$, subscript $i+\frac{1}{2}$ is the +link between node $i$ and node $i+1$, and subscript $i-\frac{3}{2}$ is the link between node +$i-1$ and node $i-2$. $\SIb{Q}{m^3 s^{-1}}$ is the water discharge, $\SIb{\eta}{m}$ is the +water surface elevation, $\SIb{h_f}{m}$ is the water depth between cells, +$\SIb{n}{m^{-\frac{1}{3}} s}$ is the Manning's roughness coefficient, $\SIb{g}{m s^{-2}}$ is +acceleration due to gravity, $\SIb{\Delta t}{s}$ is the adaptive model time step, $\SIb{\Delta +x}{m}$ is the distance between two cells and $\SIb{\Delta y}{m}$ is the flow width. Below the +staggered grid and variables of the numerical solution in the x-direction, based on Almeida et +al. (2012): ![](../../images/numerical_scheme_almeida.png) -The overland flow local inertial approach is used in combination with the local inertial -river routing. This is a similar to the modelling approach of Neal et al. (2012), where the -hydraulic model LISFLOOD-FP was extended with a subgrid channel model. For the subgrid -channel, Neal et al. (2012) make use of a D4 (four direction) scheme, while here a D8 (eight -direction) scheme is used, in combination with the D4 scheme for 2D overland flow. +The overland flow local inertial approach is used in combination with the local inertial river +routing. This is a similar to the modelling approach of Neal et al. (2012), where the hydraulic +model LISFLOOD-FP was extended with a subgrid channel model. For the subgrid channel, Neal et +al. (2012) make use of a D4 (four direction) scheme, while here a D8 (eight direction) scheme +is used, in combination with the D4 scheme for 2D overland flow. In the TOML file the following properties related to the local inertial model with 1D river routing and 2D overland flow can be provided for the `sbm` model type: @@ -143,10 +145,10 @@ froude_limit = true # default is true, limit flow to subcritical-cr h_thresh = 0.1 # water depth [m] threshold for calculating flow between cells (default = 1e-03) ``` -The properties `inertial_flow_alpha`, `froude_limit` and `h_thresh` apply to 1D river -routing as well as 2D overland flow. The properties `inertial_flow_alpha` and -`froude_limit`, and the adaptive model time step $\Delta t$ are explained in more detail -in the [River and floodplain routing](@ref) section of the local inertial model. +The properties `inertial_flow_alpha`, `froude_limit` and `h_thresh` apply to 1D river routing +as well as 2D overland flow. The properties `inertial_flow_alpha` and `froude_limit`, and the +adaptive model time step $\Delta t$ are explained in more detail in the [River and floodplain +routing](@ref) section of the local inertial model. ## Inflow External water (supply/abstraction) `inflow` $\SIb{}{m^3 s^{-1}}$ can be added to the local @@ -155,29 +157,29 @@ parameter or as part of forcing (see also [Input section](@ref)). ## Abstractions Abstractions from the river through the variable `abstraction` $\SIb{}{m^3 s^{-1}}$ are -possible when water demand and allocation is computed. The variable `abstraction` is set -from the water demand and allocation module each time step. Abstractions are subtracted as -part of the continuity equation of the local inertial model. +possible when water demand and allocation is computed. The variable `abstraction` is set from +the water demand and allocation module each time step. Abstractions are subtracted as part of +the continuity equation of the local inertial model. ## Multi-Threading -The local inertial model for river flow (1D) and river and overland flow combined (1D-2D) -can be executed in parallel using multiple threads. +The local inertial model for river flow (1D) and river and overland flow combined (1D-2D) can +be executed in parallel using multiple threads. ## References + Adams, J. M., Gasparini, N. M., Hobley, D. E. J., Tucker, G. E., Hutton, E. W. H., - Nudurupati, S. S., and Istanbulluoglu, E., 2017, The Landlab v1.0 OverlandFlow component: - a Python tool for computing shallow-water flow across watersheds, Geosci. Model Dev., 10, + Nudurupati, S. S., and Istanbulluoglu, E., 2017, The Landlab v1.0 OverlandFlow component: a + Python tool for computing shallow-water flow across watersheds, Geosci. Model Dev., 10, 1645–1663, . -+ de Almeida, G. A. M., P. Bates, J. E. Freer, and M. Souvignet, 2012, Improving the - stability of a simple formulation of the shallow water equations for 2-D flood modeling, - Water Resour. Res., 48, W05528, . -+ Bates, P. D., M. S. Horritt, and T. J. Fewtrell, 2010, A simple inertial formulation of - the shallow water equations for efficient two-dimensional flood inundation modelling, J. - Hydrol., 387, 33–45, . -+ Coulthard, T. J., Neal, J. C., Bates, P. D., Ramirez, J., de Almeida, G. A. M., and - Hancock, G. R., 2013, Integrating the LISFLOOD-FP 2- D hydrodynamic model with the CAESAR - model: implications for modelling landscape evolution, Earth Surf. Proc. Land., 38, - 1897–1906, . ++ de Almeida, G. A. M., P. Bates, J. E. Freer, and M. Souvignet, 2012, Improving the stability + of a simple formulation of the shallow water equations for 2-D flood modeling, Water Resour. + Res., 48, W05528, . ++ Bates, P. D., M. S. Horritt, and T. J. Fewtrell, 2010, A simple inertial formulation of the + shallow water equations for efficient two-dimensional flood inundation modelling, J. Hydrol., + 387, 33–45, . ++ Coulthard, T. J., Neal, J. C., Bates, P. D., Ramirez, J., de Almeida, G. A. M., and Hancock, + G. R., 2013, Integrating the LISFLOOD-FP 2- D hydrodynamic model with the CAESAR model: + implications for modelling landscape evolution, Earth Surf. Proc. Land., 38, 1897–1906, + . + Neal, J., G. Schumann, and P. Bates (2012), A subgrid channel model for simulating river - hydraulics and floodplaininundation over large and data sparse areas, Water Resour.Res., - 48, W11506, . + hydraulics and floodplaininundation over large and data sparse areas, Water Resour.Res., 48, + W11506, . diff --git a/docs/model_docs/lateral/sediment_flux.qmd b/docs/model_docs/lateral/sediment_flux.qmd index 261682394..a733b98cc 100644 --- a/docs/model_docs/lateral/sediment_flux.qmd +++ b/docs/model_docs/lateral/sediment_flux.qmd @@ -9,11 +9,11 @@ distinguished in two different structures. ## Inland Sediment Model ### Sediment Flux in overland flow -Once the amount of soil detached by both rainfall and overland flow has been estimated, it -has then to be routed and delivered to the river network. Inland routing in sediment models -is usually done by comparing the amount of detached sediment with the transport capacity of -the flow, which is the maximum amount of sediment that the flow can carry downslope. There -are several existing formulas available in the literature. For a wide range of slopes and for +Once the amount of soil detached by both rainfall and overland flow has been estimated, it has +then to be routed and delivered to the river network. Inland routing in sediment models is +usually done by comparing the amount of detached sediment with the transport capacity of the +flow, which is the maximum amount of sediment that the flow can carry downslope. There are +several existing formulas available in the literature. For a wide range of slopes and for overland flow, the Govers equation (1990) seems the most appropriate choice (Hessel et al, 2007). However, as the wflow\_sediment model was developed to be linked to water quality issues, the Yalin transport equation was chosen as it can handle particle differentiation @@ -33,7 +33,7 @@ $$ \mathrm{PSI} = 0.13\mathrm{SIL}\\ \mathrm{PCL} = 0.20\mathrm{CLA} \\ - \mathrm{SAG} = + \mathrm{SAG} = \begin{align*} \begin{cases} 2.0\mathrm{CLA} &\text{ if }\quad \mathrm{CLA} < 0.25 \\ @@ -46,25 +46,25 @@ $$ \end{gathered} $$ -where $\mathrm{CLA}$, $\mathrm{SIL}$ and $\mathrm{SAN}$ are the primary clay, silt, sand fractions of the topsoil -and $\mathrm{PCL}$, $\mathrm{PSI}$, $\mathrm{PSA}$, $\mathrm{SAG}$ and $\mathrm{LAG}$ are the clay, silt, sand, small and large -aggregates fractions of the detached sediment respectively. The transport capacity of the -flow using Yalin's equation with particle differentiation, developed by Foster (1982), is: +where $\mathrm{CLA}$, $\mathrm{SIL}$ and $\mathrm{SAN}$ are the primary clay, silt, sand +fractions of the topsoil and $\mathrm{PCL}$, $\mathrm{PSI}$, $\mathrm{PSA}$, $\mathrm{SAG}$ and +$\mathrm{LAG}$ are the clay, silt, sand, small and large aggregates fractions of the detached +sediment respectively. The transport capacity of the flow using Yalin's equation with particle +differentiation, developed by Foster (1982), is: $$ \mathrm{TC}_i = (P_e)_i (S_g)_i \, \rho_w \, g \, d_i V_* $$ where $\mathrm{TC}_i$ is the transport capacity of the flow for the particle class $i$, $(P_e)_i$ is the effective number of particles of class $i$, $\SIb{(S_g)_i}{kg m^{-3}}$ is the -specific gravity for the particle class $i$, $\SIb{\rho_w}{kg m^{-3}}$ is the mass density -of the fluid, $\SIb{g}{m s^{-2}}$ is the acceleration due to gravity, -$\SIb{d_i}{m}$ is the diameter of the particle of class $i$ and $V_* = \SIb{(g R S)^{0.5}}{m s^{-1}}$ is the -shear velocity of the flow with $S$ the slope gradient and $\SIb{R}{m}$ the -hydraulic radius of the flow. The detached sediment are then routed down slope until the -river network using the `accucapacityflux`, `accupacitystate` functions depending on the -transport capacity from Yalin. - -The choice of transport capacity method for the overland flow is set up in the model section -of the TOML: +specific gravity for the particle class $i$, $\SIb{\rho_w}{kg m^{-3}}$ is the mass density of +the fluid, $\SIb{g}{m s^{-2}}$ is the acceleration due to gravity, $\SIb{d_i}{m}$ is the +diameter of the particle of class $i$ and $V_* = \SIb{(g R S)^{0.5}}{m s^{-1}}$ is the shear +velocity of the flow with $S$ the slope gradient and $\SIb{R}{m}$ the hydraulic radius of the +flow. The detached sediment are then routed down slope until the river network using the +`accucapacityflux`, `accupacitystate` functions depending on the transport capacity from Yalin. + +The choice of transport capacity method for the overland flow is set up in the model section of +the TOML: ```toml [model] landtransportmethod = "yalinpart" # Overland flow transport capacity method: ["yalinpart", "govers", "yalin"] @@ -74,42 +74,41 @@ flow and can therefore not be used in combination with the river part of the sed ## River Sediment Model -Sediment dynamics in rivers can be described by the same three processes on -land: erosion, deposition and transport. The difference is that channel flow is much higher, -deeper and permanent compared to overland flow. In channels, erosion is the direct removal of -sediments from the river bed or bank (lateral erosion). Sediments are transported in the -river either by rolling, sliding and silting (bed load transport) or via turbulent flow in -the higher water column (suspended load transport). The type of transport is determined by -the river bed shear stress. As sediment particles have a higher density than water, they can -also be deposited on the river bed according to their settling velocity compared to the flow -velocity. In addition to regular deposition in the river, lakes, reservoirs and floodplains -represents additional major sediment settling pools. +Sediment dynamics in rivers can be described by the same three processes on land: erosion, +deposition and transport. The difference is that channel flow is much higher, deeper and +permanent compared to overland flow. In channels, erosion is the direct removal of sediments +from the river bed or bank (lateral erosion). Sediments are transported in the river either by +rolling, sliding and silting (bed load transport) or via turbulent flow in the higher water +column (suspended load transport). The type of transport is determined by the river bed shear +stress. As sediment particles have a higher density than water, they can also be deposited on +the river bed according to their settling velocity compared to the flow velocity. In addition to +regular deposition in the river, lakes, reservoirs and floodplains represents additional major +sediment settling pools. Complete models of sediment dynamics based on hydrology and not on hydraulics or hydrodynamics are much rarer than for soil loss and inland dynamics. The simpler models such as the SWAT -default sediment river model uses again the transport capacity of the flow to determine if -there is erosion or deposition (Neitsch et al., 2011). A more physics-based approach -(Partheniades, 1965) to determine river erosion is used by Liu et al. (2018) and in the new -SWAT's approach developed by Narasimhan et al. (2017). For wflow\_sediment, the new -physics-based model of SWAT was chosen for transport and erosion as it enables the use of -parameter estimation for erosion of bed and bank of the channel and separates the suspended -from the bed loads. +default sediment river model uses again the transport capacity of the flow to determine if there +is erosion or deposition (Neitsch et al., 2011). A more physics-based approach (Partheniades, +1965) to determine river erosion is used by Liu et al. (2018) and in the new SWAT's approach +developed by Narasimhan et al. (2017). For wflow\_sediment, the new physics-based model of SWAT +was chosen for transport and erosion as it enables the use of parameter estimation for erosion +of bed and bank of the channel and separates the suspended from the bed loads. ![Overview of the different processes for a river cell in wflow\_sediment.](../../images/river-scheme.png) -Running the river model is an option of the wflow\_sediment model and is enabled using the -TOML file. By default it is `false`: +Running the river model is an option of the wflow\_sediment model and is enabled using the TOML +file. By default it is `false`: ```toml [model] runrivermodel = true ``` ### Sediment inputs in a river cell -The first part of the river model assesses how much detached sediment are in the river cell -at the beginning of the timestep $t$. Sources of detached sediment are sediments coming -from land erosion, estimated with the soil loss part of wflow_sediment model, the sediment -coming from upstream river cells and the detached sediment that were left in the cell at the -end of the previous timestep ``(t-1)``: +The first part of the river model assesses how much detached sediment are in the river cell at +the beginning of the timestep $t$. Sources of detached sediment are sediments coming from land +erosion, estimated with the soil loss part of wflow_sediment model, the sediment coming from +upstream river cells and the detached sediment that were left in the cell at the end of the +previous timestep ``(t-1)``: $$ (\subtext{\mathrm{sed}}{in})_t = (\subtext{\mathrm{sed}}{land})_t + \mathrm{upstream}\left[(\subtext{\text{sed}}{out})_{t-1}\right] + (\subtext{\text{sed}}{riv})_{t-1} $$ @@ -119,11 +118,10 @@ Once the amount of sediment inputs at the beginning of the timestep is known, th estimates transport, and river erosion if there is a deficit of sediments. Transport in the river system is estimated via a transport capacity formula. There are several transport capacity formulas available in wflow_sediment, some requiring calibration and some not. -Choosing a transport capacity equation depends on the river characteristics (some equation -are more suited for narrow or wider rivers), and on the reliability of the required river -parameters (such as slope, width or mean particle diameter of the river channel). Several -river transport capacity are available and the choice is set up in the model section of the -TOML: +Choosing a transport capacity equation depends on the river characteristics (some equation are +more suited for narrow or wider rivers), and on the reliability of the required river +parameters (such as slope, width or mean particle diameter of the river channel). Several river +transport capacity are available and the choice is set up in the model section of the TOML: ```toml [model] rivtransportmethod = "bagnold" # River flow transport capacity method: ["bagnold", "engelund", "yang", "kodatie", "molinas"] @@ -131,21 +129,23 @@ rivtransportmethod = "bagnold" # River flow transport capacity method: ["bagnold **Simplified Bagnold** -Originally more valid for intermediate to large rivers, this simplified version of the -Bagnold equation relates sediment transport to flow velocity with two simple calibration -parameters (Neitsch et al, 2011): +Originally more valid for intermediate to large rivers, this simplified version of the Bagnold +equation relates sediment transport to flow velocity with two simple calibration parameters +(Neitsch et al, 2011): $$ C_{\max} = \subtext{c}{sp} \left( \dfrac{\mathrm{prf} Q}{h W} \right)^{\subtext{\mathrm{sp}}{exp}} $$ -where $\SIb{C_{\max}}{kg L^{-1}}$ (or $\SIb{}{ton m^{-1}}$) is the sediment concentration, $\SIb{Q}{m^3 s^{-1}}$ is the -surface runoff in the river cell, $\SIb{h}{m}$ is the river water level, -$\SIb{W}{m}$ is the river width and $\subtext{c}{sp}$, $\mathrm{prf}$ and $\subtext{\mathrm{sp}}{exp}$ are calibration -parameters. The $\mathrm{prf}$ coefficient is usually used to deduce the peak velocity of the flow, -but for simplification in `wflow_sediment`, the equation was simplified to only get two -parameters to calibrate: $\subtext{\mathrm{sp}}{exp}$ and $\subtext{c}{Bagnold} = \subtext{c}{sp} \, \mathrm{prf}^{\subtext{\mathrm{sp}}{exp}}$. The -coefficient $\subtext{\mathrm{sp}}{exp}$ usually varies between $1$ and $2$ while $\mathrm{prf}$ and $\subtext{c}{sp}$ have a -wider range of variation. The table below summarizes ranges and values of the three Bagnold +where $\SIb{C_{\max}}{kg L^{-1}}$ (or $\SIb{}{ton m^{-1}}$) is the sediment concentration, +$\SIb{Q}{m^3 s^{-1}}$ is the surface runoff in the river cell, $\SIb{h}{m}$ is the river water +level, $\SIb{W}{m}$ is the river width and $\subtext{c}{sp}$, $\mathrm{prf}$ and +$\subtext{\mathrm{sp}}{exp}$ are calibration parameters. The $\mathrm{prf}$ coefficient is +usually used to deduce the peak velocity of the flow, but for simplification in +`wflow_sediment`, the equation was simplified to only get two parameters to calibrate: +$\subtext{\mathrm{sp}}{exp}$ and $\subtext{c}{Bagnold} = \subtext{c}{sp} \, +\mathrm{prf}^{\subtext{\mathrm{sp}}{exp}}$. The coefficient $\subtext{\mathrm{sp}}{exp}$ +usually varies between $1$ and $2$ while $\mathrm{prf}$ and $\subtext{c}{sp}$ have a wider +range of variation. The table below summarizes ranges and values of the three Bagnold coefficients used by other studies: Table: Range of the simplified Bagnold coefficients (and calibrated value) @@ -157,16 +157,17 @@ Table: Range of the simplified Bagnold coefficients (and calibrated value) | Abbaspour 2007 | Thur (CH) | 0.2-0.25 (/) | 0.001-0.002 (/) | 0.35-1.47 (/) | | Oeurng 2011 | Save (FR) | 0-2 (0.58) | 0.0001-0.01 (0.01) | 1-2 (2) | -**Engelund and Hansen** This transport capacity is not present in SWAT but used in many -models such as Delft3D-WAQ, Engelund and Hansen calculates the total sediment load as -(Engelund and Hansen, 1967): +**Engelund and Hansen** This transport capacity is not present in SWAT but used in many models +such as Delft3D-WAQ, Engelund and Hansen calculates the total sediment load as (Engelund and +Hansen, 1967): $$ C_w = 0.05 \left( \dfrac{\rho_{s}}{\rho_{s} - \rho} \right) \left( \dfrac{u S}{\sqrt{\left( \dfrac{\rho_{s}}{\rho_{s} - \rho} \right) g D_{50}}} \right) \theta^{1/2} $$ where $C_w$ is the sediment concentration by weight, $\SIb{\rho}{g m^{-3}}$ and $\SIb{\rho_{s}}{g m^{-3}}$ are the -fluid and sediment density (here respectively equal to $\SI{1000}{g m^{-3}}$ and $\SI{2650}{g m^{-3}}$), $\SIb{u}{m s^{-1}}$ is the water -mean velocity, $S$ is the river slope, $g$ is the acceleration due to gravity, -$\SIb{D_{50}}{m}$ is the river mean diameter and $\theta$ is the Shields parameter. +fluid and sediment density (here respectively equal to $\SI{1000}{g m^{-3}}$ and $\SI{2650}{g +m^{-3}}$), $\SIb{u}{m s^{-1}}$ is the water mean velocity, $S$ is the river slope, $g$ is the +acceleration due to gravity, $\SIb{D_{50}}{m}$ is the river mean diameter and $\theta$ is the +Shields parameter. **Kodatie** Kodatie (1999) developed the power relationships from Posada (1995) using field data and linear optimization so that they would be applicable for a wider range of riverbed @@ -175,9 +176,8 @@ $$ C_{\max} = \left( \dfrac{a u^{b} h^{c} S^{d}}{\subtext{V}{in}} \right) W $$ where $\SIb{\subtext{V}{in}}{m^3}$ in the volume of water entering the river cell -during the timestep and $a$, $b$, $c$ and $d$ are coefficients depending -on the riverbed sediment size. Values of these coefficients are summarized in the table -below. +during the timestep and $a$, $b$, $c$ and $d$ are coefficients depending on the riverbed +sediment size. Values of these coefficients are summarized in the table below. Table: Range of the simplified Bagnold coefficients (and calibrated value) @@ -200,9 +200,10 @@ $$ +\left(2.784-0.305\log\left(\frac{\omega_{s,50}D_{50}}{\nu}\right)-0.282\log\left(\frac{u_*}{\omega_{s,50}}\right)\right)\log\left(\frac{uS}{\omega_{s,50}}-\frac{u_{cr}S}{\omega_{s,50}}\right) $$ where $C_{ppm}$ is sediment concentration in parts per million by weight, -$\SIb{\omega_{s,50}}{m s^{-1}}$ is the settling velocity of a particle with the median riverbed diameter -estimated with Stokes, $\SIb{\nu}{m^2 s^{-1}}$ is the kinematic viscosity of the fluid, $\SIb{u_*}{m s^{-1}}$ is the shear velocity where $u_* = \sqrt{gR_{H}S}$ -with $R_{H}$ the hydraulic radius of the river and $\SIb{u_{cr}}{m s^{-1}}$ is the critical velocity +$\SIb{\omega_{s,50}}{m s^{-1}}$ is the settling velocity of a particle with the median riverbed +diameter estimated with Stokes, $\SIb{\nu}{m^2 s^{-1}}$ is the kinematic viscosity of the +fluid, $\SIb{u_*}{m s^{-1}}$ is the shear velocity where $u_* = \sqrt{gR_{H}S}$ with $R_{H}$ +the hydraulic radius of the river and $\SIb{u_{cr}}{m s^{-1}}$ is the critical velocity (equation can be found in Hessel, 2007). **Molinas and Wu** The Molinas and Wu (2001) transport equation was developed for large @@ -217,25 +218,28 @@ $$ $$ Once the maximum concentration $C_{\max}$ is established with one of the above transport -formula, the model then determines if there is erosion of the river bed and bank. In order -to do that, the difference $sed_{ex}$ between the maximum amount of sediment estimated -with transport ($\mathrm{sed}_{\max} = C_{\max} \subtext{V}{in}$) and the sediment inputs to the river cell -($\subtext{\mathrm{sed}}{in}$ calculated above) is calculated. If too much sediment is coming in and -$\subtext{\mathrm{sed}}{ex}$ is negative, then there is no river bed and bank erosion. And if the river has -not reach its maximum transport capacity, then erosion of the river happens. - -First, the sediments stored in the cell from deposition in previous timesteps $\subtext{\mathrm{sed}}{stor}$ -are eroded from clay to gravel. If this amount is not enough to cover $\subtext{\mathrm{sed}}{ex}$, then -erosion of the local river bed and bank material starts. +formula, the model then determines if there is erosion of the river bed and bank. In order to +do that, the difference $sed_{ex}$ between the maximum amount of sediment estimated with +transport ($\mathrm{sed}_{\max} = C_{\max} \subtext{V}{in}$) and the sediment inputs to the +river cell ($\subtext{\mathrm{sed}}{in}$ calculated above) is calculated. If too much sediment +is coming in and $\subtext{\mathrm{sed}}{ex}$ is negative, then there is no river bed and bank +erosion. And if the river has not reach its maximum transport capacity, then erosion of the +river happens. + +First, the sediments stored in the cell from deposition in previous timesteps +$\subtext{\mathrm{sed}}{stor}$ are eroded from clay to gravel. If this amount is not enough to +cover $\subtext{\mathrm{sed}}{ex}$, then erosion of the local river bed and bank material +starts. Instead of just setting river erosion amount to just cover the remaining difference -$\subtext{\mathrm{sed}}{exeff}$ between $\subtext{\mathrm{sed}}{ex}$ and $\subtext{\mathrm{sed}}{stor}$, actual erosion potential is -adjusted using river characteristics and is separated between the bed and bank of the river -using the physics-based approach of Knight (1984). - -The bed and bank of the river are supposed to only be able to erode a maximum amount of -their material ``E_{R,\mathrm{bed}}`` for the bed and ``E_{R,\mathrm{bank}}`` for the river bank. For a -rectangular channel, assuming it is meandering and thus only one bank is prone to erosion, +$\subtext{\mathrm{sed}}{exeff}$ between $\subtext{\mathrm{sed}}{ex}$ and +$\subtext{\mathrm{sed}}{stor}$, actual erosion potential is adjusted using river +characteristics and is separated between the bed and bank of the river using the physics-based +approach of Knight (1984). + +The bed and bank of the river are supposed to only be able to erode a maximum amount of their +material ``E_{R,\mathrm{bed}}`` for the bed and ``E_{R,\mathrm{bank}}`` for the river bank. For +a rectangular channel, assuming it is meandering and thus only one bank is prone to erosion, they are calculated from the equations (Neitsch et al, 2011): $$ \begin{gathered} @@ -243,39 +247,39 @@ $$ E_{R,\mathrm{bank}} = k_{d,\mathrm{bank}} \left( \tau_{e,\mathrm{bank}} - \tau_{cr,\mathrm{bank}} \right) 10^{-6} L h \rho_{b, \mathrm{bank}} \Delta t \end{gathered} $$ -where $\SIb{E_R}{ton}$ is the potential bed/bank erosion rates, $\SIb{k_d}{cm^3 N^{-1}, s^{-1}}$ is the erodibility -of the bed/bank material, $\SIb{\tau_e}{N m^{-2}}$ is the effective -shear stress from the flow on the bed/bank, $\SIb{\tau_{cr}}{N m^{-2}}$ is the critical -shear stress for erosion to happen, $\SIb{L}{m}$, $\SIb{W}{m}$ and $\SIb{h}{m}$ are the channel -length, width and water height, $\SIb{\rho_{b}}{g cm^{-3}}$ is the bulk density of the bed/bank of the -river and $\SIb{\Delta t}{s}$ is the model timestep. +where $\SIb{E_R}{ton}$ is the potential bed/bank erosion rates, $\SIb{k_d}{cm^3 N^{-1}, +s^{-1}}$ is the erodibility of the bed/bank material, $\SIb{\tau_e}{N m^{-2}}$ is the effective +shear stress from the flow on the bed/bank, $\SIb{\tau_{cr}}{N m^{-2}}$ is the critical shear +stress for erosion to happen, $\SIb{L}{m}$, $\SIb{W}{m}$ and $\SIb{h}{m}$ are the channel +length, width and water height, $\SIb{\rho_{b}}{g cm^{-3}}$ is the bulk density of the bed/bank +of the river and $\SIb{\Delta t}{s}$ is the model timestep. In wflow_sediment, the erodibility of the bed/bank are approximated using the formula from Hanson and Simon (2001): $$ k_d=0.2 \tau_{cr}^{-0.5} $$ -Normally erodibilities are evaluated using jet test in the field and there are several -reviews and some adjustments possible to this equation (Simon et al, 2011). However, to -avoid too heavy calibration and for the scale considered, this equation is supposed to be -efficient enough. The critical shear stress $\tau_{cr}$ is evaluated differently for the -bed and bank. For the bed, the most common formula from Shields initiation of movement is -used. For the bank, a more recent approach from Julian and Torres (2006) is used : +Normally erodibilities are evaluated using jet test in the field and there are several reviews +and some adjustments possible to this equation (Simon et al, 2011). However, to avoid too heavy +calibration and for the scale considered, this equation is supposed to be efficient enough. The +critical shear stress $\tau_{cr}$ is evaluated differently for the bed and bank. For the bed, +the most common formula from Shields initiation of movement is used. For the bank, a more +recent approach from Julian and Torres (2006) is used : $$ \tau_{cr,\mathrm{bank}} = (0.1+0.1779 SC+0.0028 SC^{2}-2.34 10^{-5} SC^{3}) C_{ch} $$ -where $SC$ is the percent clay and silt content of the river bank and $C_{ch}$ is a -coefficient taking into account the positive impact of vegetation on erosion reduction. This -coefficient is then dependent on the land use and classical values are shown in the table -below. These values where then adapted for use with the GlobCover land use map. Percent of -clay and silt (along with sand and gravel) for the channel is estimated from the river -median particle diameter assuming the same values as SWAT shown in the table below. Median -particle diameter is here estimated depending on the Strahler river order. The higher the -order, the smaller the diameter is. As the median diameter is only used in wflow_sediment -for the estimation of the river bed/bank sediment composition, this supposition should be -enough. Actual refined data or calibration may however be needed if the median diameter is -also required for the transport formula. In a similar way, the bulk densities of river bed -and bank are also just assumed to be of respectively $\SI{1.5}{g cm^{-3}}$ and $\SI{1.4}{g cm^{-3}}$. +where $SC$ is the percent clay and silt content of the river bank and $C_{ch}$ is a coefficient +taking into account the positive impact of vegetation on erosion reduction. This coefficient is +then dependent on the land use and classical values are shown in the table below. These values +where then adapted for use with the GlobCover land use map. Percent of clay and silt (along +with sand and gravel) for the channel is estimated from the river median particle diameter +assuming the same values as SWAT shown in the table below. Median particle diameter is here +estimated depending on the Strahler river order. The higher the order, the smaller the diameter +is. As the median diameter is only used in wflow_sediment for the estimation of the river +bed/bank sediment composition, this supposition should be enough. Actual refined data or +calibration may however be needed if the median diameter is also required for the transport +formula. In a similar way, the bulk densities of river bed and bank are also just assumed to be +of respectively $\SI{1.5}{g cm^{-3}}$ and $\SI{1.4}{g cm^{-3}}$. Table: Classical values of the channel cover vegetation coefficient (Julian and Torres, 2006) @@ -286,7 +290,8 @@ Table: Classical values of the channel cover vegetation coefficient (Julian and | Sparse trees | 5.40 | | Dense trees | 19.20 | -Table: Composition of the river bed/bank depending on the median diameter $\SIb{d_{50}}{\mu m}$ (Neitsch et al, 2011) +Table: Composition of the river bed/bank depending on the median diameter $\SIb{d_{50}}{\mu m}$ +(Neitsch et al, 2011) |Sediment Fraction | $\leq$ 5 | 5 to 50 | 50 to 2000 | $>$ 2000 | | ---- | ---- | ---- | ---- | ---- | @@ -295,18 +300,19 @@ Table: Composition of the river bed/bank depending on the median diameter $\SIb{ | Clay | 0.65 | 0.15 | 0.15 | 0.05 | | Gravel | 0.05 | 0.05 | 0.05 | 0.65 | -Then, the repartition of the flow shear stress is refined into the effective shear stress -and the bed and bank of the river using the equations developed by Knight (1984) for a -rectangular channel: +Then, the repartition of the flow shear stress is refined into the effective shear stress and +the bed and bank of the river using the equations developed by Knight (1984) for a rectangular +channel: $$ \begin{gathered} \tau_{e,\mathrm{bed}} = \rho g R_{H} S \left(1 - \dfrac{SF_{\mathrm{bank}}}{100}\right) \left(1+\dfrac{2h}{W}\right) \\ \tau_{e,\mathrm{bank}} = \rho g R_{H} S \left( SF_{\mathrm{bank}}\right) \left(1+\dfrac{W}{2h}\right) \end{gathered} $$ -where $\rho g$ is the fluid specific weight ($\SI{9800}{N m^{-3}}$ for water), $\SIb{R_H}{m}$ is the -hydraulic radius of the channel, $\SIb{h}{m}$ and $\SIb{W}{m}$ are the water level and river width. $SF_{\mathrm{bank}}$ is the proportion of shear stress acting on the bank (%) and is estimated -from (Knight, 1984): +where $\rho g$ is the fluid specific weight ($\SI{9800}{N m^{-3}}$ for water), $\SIb{R_H}{m}$ +is the hydraulic radius of the channel, $\SIb{h}{m}$ and $\SIb{W}{m}$ are the water level and +river width. $SF_{\mathrm{bank}}$ is the proportion of shear stress acting on the bank (%) and +is estimated from (Knight, 1984): $$ \mathrm{SF}_{\mathrm{bank}} = \exp \left( -3.230 \log_{10}\left(\dfrac{W}{h}+3\right)+6.146 \right) $$ @@ -319,8 +325,8 @@ $$ $$ And the final actual eroded amount for the bed and bank is the maximum between $\mathrm{RTE} \subtext{\mathrm{sed}}{exeff}$ and the erosion potential $E_R$. Total eroded amount of sediment -$\subtext{\mathrm{sed}}{erod}$ is then the sum of the eroded sediment coming from the storage of previously -deposited sediment and the river bed/bank erosion. +$\subtext{\mathrm{sed}}{erod}$ is then the sum of the eroded sediment coming from the storage +of previously deposited sediment and the river bed/bank erosion. ### River deposition As sediments have a higher density than water, moving sediments in water can be deposited in @@ -330,17 +336,18 @@ Einstein's equation (Neitsch et al, 2011): $$ \subtext{P}{dep}=\left(1-\dfrac{1}{e^{x}}\right)100 $$ -where $\subtext{P}{dep}$ is the percentage of sediments that is deposited on the river bed and x is -a parameter calculated with: +where $\subtext{P}{dep}$ is the percentage of sediments that is deposited on the river bed and +x is a parameter calculated with: $$ x = \dfrac{1.055 L \omega_{s}}{u h} $$ -where $\SIb{L}{m}$ and $\SIb{h}{m}$ are channel length and water height, $\SIb{\omega_s}{m s^{-1}}$ is the -particle settling velocity calculated with Stokes' formula and $\SIb{u}{m s^{-1}}$ is the -mean flow velocity. The calculated percentage is then subtracted from the -amount of sediment input and eroded river sediment for each particle size class ($\subtext{\mathrm{sed}}{dep} -= \subtext{P}{dep}/100 (\subtext{\mathrm{sed}}{in} + \subtext{\mathrm{sed}}{erod})$). Resulting deposited sediment are then stored in -the river bed and can be re-mobilized in future time steps by erosion. +where $\SIb{L}{m}$ and $\SIb{h}{m}$ are channel length and water height, $\SIb{\omega_s}{m +s^{-1}}$ is the particle settling velocity calculated with Stokes' formula and $\SIb{u}{m +s^{-1}}$ is the mean flow velocity. The calculated percentage is then subtracted from the +amount of sediment input and eroded river sediment for each particle size class +($\subtext{\mathrm{sed}}{dep} = \subtext{P}{dep}/100 (\subtext{\mathrm{sed}}{in} + +\subtext{\mathrm{sed}}{erod})$). Resulting deposited sediment are then stored in the river bed +and can be re-mobilized in future time steps by erosion. ### Mass balance and sediment concentration Finally after estimating inputs, deposition and erosion with the transport capacity of the @@ -350,13 +357,14 @@ $$ \subtext{\mathrm{sed}}{out} = (\subtext{\mathrm{sed}}{in} + \subtext{\mathrm{sed}}{erod} - \subtext{\mathrm{sed}}{dep}) \dfrac{\subtext{V}{out}}{V} $$ -where $\SIb{\subtext{\mathrm{sed}}{out}}{ton}$ is the amount of sediment leaving the river cell (tons), $\SIb{\subtext{\mathrm{sed}}{in}}{ton}$ is -the amount of sediment coming into the river cell (storage from previous timestep, land -erosion and sediment flux from upstream river cells), $\SIb{\subtext{\mathrm{sed}}{erod}}{ton}$ is the amount -of sediment coming from river erosion, $\SIb{\subtext{\mathrm{sed}}{dep}}{ton}$ is the amount of deposited -sediments, $\SIb{\subtext{V}{out}}{m^3}$ is the volume of water leaving the river cell (surface runoff -$Q$ times timestep $\Delta t$) and $\SIb{V}{m^3}$ is the total volume of water in -the river cell ($\subtext{V}{out}$ plus storage $h W L$). +where $\SIb{\subtext{\mathrm{sed}}{out}}{ton}$ is the amount of sediment leaving the river cell +(tons), $\SIb{\subtext{\mathrm{sed}}{in}}{ton}$ is the amount of sediment coming into the river +cell (storage from previous timestep, land erosion and sediment flux from upstream river +cells), $\SIb{\subtext{\mathrm{sed}}{erod}}{ton}$ is the amount of sediment coming from river +erosion, $\SIb{\subtext{\mathrm{sed}}{dep}}{ton}$ is the amount of deposited sediments, +$\SIb{\subtext{V}{out}}{m^3}$ is the volume of water leaving the river cell (surface runoff $Q$ +times timestep $\Delta t$) and $\SIb{V}{m^3}$ is the total volume of water in the river cell +($\subtext{V}{out}$ plus storage $h W L$). A mass balance is then used to calculate the amount of sediment remaining in the cell at the end of the timestep $(\subtext{\mathrm{sed}}{riv})_t$: @@ -367,30 +375,31 @@ $$ ### Lake and reservoir modelling Apart from land and river, the hydrologic `wflow_sbm` model also handles lakes and reservoirs modelling. In `wflow_sbm`, lakes and large reservoirs are modelled using a 1D bucket model at -the cell corresponding to the outlet. For the other cells belonging to the lake/reservoir -which are not the outlet, processes such as precipitation and evaporation are filtered out -and shifted to the outlet cell. `wflow_sediment` handles the lakes and reservoirs in the same way. If a -cell belongs to a lake/reservoir and is not the outlet then the model assumes that no -erosion/deposition of sediments is happening and the sediments are only all transported to -the lake/reservoir outlet. Once the sediments reach the outlet, then sediments are deposited -in the lake/reservoir according to Camp's model (1945) (Verstraeten et al, 2000): +the cell corresponding to the outlet. For the other cells belonging to the lake/reservoir which +are not the outlet, processes such as precipitation and evaporation are filtered out and +shifted to the outlet cell. `wflow_sediment` handles the lakes and reservoirs in the same way. +If a cell belongs to a lake/reservoir and is not the outlet then the model assumes that no +erosion/deposition of sediments is happening and the sediments are only all transported to the +lake/reservoir outlet. Once the sediments reach the outlet, then sediments are deposited in the +lake/reservoir according to Camp's model (1945) (Verstraeten et al, 2000): $$ \mathrm{TE} = \dfrac{\omega_s}{u_{cr,\mathrm{res}}} = \dfrac{\subtext{A}{res}}{\subtext{Q}{out,res}} \omega_s $$ -where $\mathrm{TE}$ is the trapping efficiency of the lake/reservoir (or the fraction of particles -trapped), $\SIb{\omega_{s}}{m s^{-1}}$ is the particle velocity from Stokes, $\SIb{\subtext{u}{cr,res}}{m s^{-1}}$ -is the reservoir's critical settling velocity which is equal to the reservoir's -outflow $\SIb{\subtext{Q}{out,res}}{m^3 s^{-1}}$ divided by the reservoir's surface area -$\SIb{\subtext{A}{res}}{m^2}$. +where $\mathrm{TE}$ is the trapping efficiency of the lake/reservoir (or the fraction of +particles trapped), $\SIb{\omega_{s}}{m s^{-1}}$ is the particle velocity from Stokes, +$\SIb{\subtext{u}{cr,res}}{m s^{-1}}$ is the reservoir's critical settling velocity which is +equal to the reservoir's outflow $\SIb{\subtext{Q}{out,res}}{m^3 s^{-1}}$ divided by the +reservoir's surface area $\SIb{\subtext{A}{res}}{m^2}$. -For reservoirs, coarse sediment particles from the bed load are also assumed to be trapped by the -dam structure. This adding trapping is taken into account with a reservoir trapping efficiency coefficient -for large particles (between $0$ and $1$). Depending on the type of the dam, all bed load particles are trapped -(`restrapefficiency = 1.0`, for example for a gravity dam) or only partly (for example for run-of-the-river dams). +For reservoirs, coarse sediment particles from the bed load are also assumed to be trapped by +the dam structure. This adding trapping is taken into account with a reservoir trapping +efficiency coefficient for large particles (between $0$ and $1$). Depending on the type of the +dam, all bed load particles are trapped (`restrapefficiency = 1.0`, for example for a gravity +dam) or only partly (for example for run-of-the-river dams). -Lake and reservoir modelling is enabled in the model section of the TOML and require the -extra following input arguments: +Lake and reservoir modelling is enabled in the model section of the TOML and require the extra +following input arguments: ```toml [model] @@ -420,14 +429,13 @@ transport in overland flow. ## References + K.C. Abbaspour, J. Yang, I. Maximov, R. Siber, K. Bogner, J. Mieleitner, J. Zobrist, and - R.Srinivasan. Modelling hydrology and water quality in the pre-alpine/alpine Thur - watershed using SWAT. Journal of Hydrology, 333(2-4):413-430, 2007. - 10.1016/j.jhydrol.2006.09.014 -+ P. Borrelli, M. Märker, P. Panagos, and B. Schütt. Modeling soil erosion and river - sediment yield for an intermountain drainage basin of the Central Apennines, Italy. - Catena, 114:45-58, 2014. 10.1016/j.catena.2013.10.007 -+ F. Engelund and E. Hansen. A monograph on sediment transport in alluvial streams. - Technical University of Denmark 0stervoldgade 10, Copenhagen K., 1967. + R.Srinivasan. Modelling hydrology and water quality in the pre-alpine/alpine Thur watershed + using SWAT. Journal of Hydrology, 333(2-4):413-430, 2007. 10.1016/j.jhydrol.2006.09.014 ++ P. Borrelli, M. Märker, P. Panagos, and B. Schütt. Modeling soil erosion and river sediment + yield for an intermountain drainage basin of the Central Apennines, Italy. Catena, 114:45-58, + 2014. 10.1016/j.catena.2013.10.007 ++ F. Engelund and E. Hansen. A monograph on sediment transport in alluvial streams. Technical + University of Denmark 0stervoldgade 10, Copenhagen K., 1967. + G. Govers. Empirical relationships for the transport capacity of overland flow. IAHS Publication, (January 1990):45-63 ST, 1990. + G.J Hanson and A Simon. Erodibility of cohesive streambeds in the loess area of the @@ -435,25 +443,24 @@ transport in overland flow. + R Hessel and V Jetten. Suitability of transport equations in modelling soil erosion for a small Loess Plateau catchment. Engineering Geology, 91(1):56-71, 2007. 10.1016/j.enggeo.2006.12.013 -+ J.P Julian, and R. Torres. Hydraulic erosion of cohesive riverbanks. Geomorphology, - 76:193-206, 2006. 10.1016/j.geomorph.2005.11.003 -+ D.W. Knight, J.D. Demetriou, and M.E. Hamed. Boundary Shear in Smooth Rectangular - Channels. J. Hydraul. Eng., 110(4):405-422, 1984. 10.1061/(ASCE)0733-9429(1987)113:1(120) -+ S.L Neitsch, J.G Arnold, J.R Kiniry, and J.R Williams. SWAT Theoretical Documentation - Version 2009. Texas Water Resources Institute, pages 1-647, 2011. - 10.1016/j.scitotenv.2015.11.063 ++ J.P Julian, and R. Torres. Hydraulic erosion of cohesive riverbanks. Geomorphology, + 76:193-206, 2006. 10.1016/j.geomorph.2005.11.003 ++ D.W. Knight, J.D. Demetriou, and M.E. Hamed. Boundary Shear in Smooth Rectangular Channels. + J. Hydraul. Eng., 110(4):405-422, 1984. 10.1061/(ASCE)0733-9429(1987)113:1(120) ++ S.L Neitsch, J.G Arnold, J.R Kiniry, and J.R Williams. SWAT Theoretical Documentation Version + 2009. Texas Water Resources Institute, pages 1-647, 2011. 10.1016/j.scitotenv.2015.11.063 + C. Oeurng, S. Sauvage, and J.M. Sanchez-Perez. Assessment of hydrology, sediment and particulate organic carbon yield in a large agricultural catchment using the SWAT model. Journal of Hydrology, 401:145-153, 2011. 10.1016/j.hydrol.2011.02.017 -+ A. Simon, N. Pollen-Bankhead, and R.E Thomas. Development and application of a - deterministic bank stability and toe erosion model for stream restoration. Geophysical - Monograph Series, 194:453-474, 2011. 10.1029/2010GM001006 ++ A. Simon, N. Pollen-Bankhead, and R.E Thomas. Development and application of a deterministic + bank stability and toe erosion model for stream restoration. Geophysical Monograph Series, + 194:453-474, 2011. 10.1029/2010GM001006 + G. Verstraeten and J. Poesen. Estimating trap efficiency of small reservoirs and ponds: methods and implications for the assessment of sediment yield. Progress in Physical Geography, 24(2):219-251, 2000. 10.1177/030913330002400204 + O. Vigiak, A. Malago, F. Bouraoui, M. Vanmaercke, and J. Poesen. Adapting SWAT hillslope - erosion model to predict sediment concentrations and yields in large Basins. Science of - the Total Environment, 538:855-875, 2015. 10.1016/j.scitotenv.2015.08.095 + erosion model to predict sediment concentrations and yields in large Basins. Science of the + Total Environment, 538:855-875, 2015. 10.1016/j.scitotenv.2015.08.095 + O. Vigiak, A. Malago, F. Bouraoui, M. Vanmaercke, F. Obreja, J. Poesen, H. Habersack, J. - Feher, and S. Groselj. Modelling sediment fluxes in the Danube River Basin with SWAT. - Science of the Total Environment, 2017. 10.1016/j.scitotenv.2017.04.236 + Feher, and S. Groselj. Modelling sediment fluxes in the Danube River Basin with SWAT. Science + of the Total Environment, 2017. 10.1016/j.scitotenv.2017.04.236 diff --git a/docs/model_docs/lateral/waterbodies.qmd b/docs/model_docs/lateral/waterbodies.qmd index d218018c8..efdf0dd91 100644 --- a/docs/model_docs/lateral/waterbodies.qmd +++ b/docs/model_docs/lateral/waterbodies.qmd @@ -5,8 +5,8 @@ title: Reservoirs and Lakes Simplified reservoirs and lakes models can be included as part of the river network. ## Reservoirs -Simple reservoirs can be included within the river routing by supplying the following -reservoir parameters: +Simple reservoirs can be included within the river routing by supplying the following reservoir +parameters: + `locs` - Outlet of the reservoirs in which each reservoir has a unique id + `area` - Surface area of the reservoirs $\SIb{}{m^2}$ @@ -25,8 +25,8 @@ lines in the TOML file of the model: [model] reservoirs = true ``` -Finally there is a mapping required between external and internal parameter names in the -TOML file, with below an example: +Finally there is a mapping required between external and internal parameter names in the TOML +file, with below an example: ```toml [input] @@ -48,28 +48,31 @@ $$ \dfrac{S(t + \Delta t)}{\Delta t} = \dfrac{S(t)}{\Delta t} + \subtext{Q}{in} + \dfrac{(P-E) A}{\Delta t} - \subtext{Q}{out} $$ -where $\SIb{S}{m^3}$ is lake storage, $\SIb{\Delta t}{s}$ is the model timestep, $\SIb{\subtext{Q}{in}}{m^3 s^{-1}}$ is -the sum of inflows (river, overland and lateral subsurface flow), -$\SIb{\subtext{Q}{out}}{m^3 s^{-1}}$ is the lake outflow at the outlet, $\SIb{P}{m}$ is precipitation, $\SIb{E}{m}$ is lake evaporation and $\SIb{A}{m^2}$ is the lake surface area. +where $\SIb{S}{m^3}$ is lake storage, $\SIb{\Delta t}{s}$ is the model timestep, +$\SIb{\subtext{Q}{in}}{m^3 s^{-1}}$ is the sum of inflows (river, overland and lateral +subsurface flow), $\SIb{\subtext{Q}{out}}{m^3 s^{-1}}$ is the lake outflow at the outlet, +$\SIb{P}{m}$ is precipitation, $\SIb{E}{m}$ is lake evaporation and $\SIb{A}{m^2}$ is the lake +surface area. ![Lake schematization.](../../images/lake.png) Most of the variables in this equation are already known or coming from previous timestep, apart from $S(t+ \Delta t)$ and $\subtext{Q}{out}$ which can both be linked to the water level -$H$ in the lake using a storage curve $S = f(H)$ and a rating curve $Q = f(H)$. In -wflow, several options are available to select storage and rating curves, and in most cases, -the mass balance is then solved by linearization and iteration or using the Modified Puls -Approach from Maniak (Burek et al., 2013). Storage curves in wflow can either: +$H$ in the lake using a storage curve $S = f(H)$ and a rating curve $Q = f(H)$. In wflow, +several options are available to select storage and rating curves, and in most cases, the mass +balance is then solved by linearization and iteration or using the Modified Puls Approach from +Maniak (Burek et al., 2013). Storage curves in wflow can either: + Come from the interpolation of field data linking volume and lake height, + Be computed from the simple relationship $S = A H$. Rating curves in wflow can either: -+ Come from the interpolation of field data linking lake outflow and water height, also appropriate for regulated lakes/ dams, -+ Be computed from a rating curve of the form $\subtext{Q}{out} = \alpha (H-H_0)^\beta$, - where $H_{0}$ is the minimum water level under which the outflow is zero. Usual values - for $\beta$ are $\frac{3}{2}$ for a rectangular weir or $2$ for a parabolic weir (Bos, 1989). ++ Come from the interpolation of field data linking lake outflow and water height, also + appropriate for regulated lakes/ dams, ++ Be computed from a rating curve of the form $\subtext{Q}{out} = \alpha (H-H_0)^\beta$, where + $H_{0}$ is the minimum water level under which the outflow is zero. Usual values for $\beta$ + are $\frac{3}{2}$ for a rectangular weir or $2$ for a parabolic weir (Bos, 1989). ### Modified Puls Approach The Modified Puls Approach is a resolution method of the lake balance that uses an explicit @@ -90,7 +93,7 @@ $$ The solution for $Q$ is then: $$ - Q = + Q = \begin{cases} \begin{align*} \frac{1}{4}\left(-\mathrm{LF} + \sqrt{\mathrm{LF}^{2} + 4 \left(\mathrm{SI} - \dfrac{A H_0}{\Delta t} \right)} @@ -116,16 +119,16 @@ following parameters: + `linkedlakelocs` - Outlet of linked (downstream) lakes (unique id) + `waterlevel` - Lake water level [m], used to reinitiate lake model + `threshold` - Water level threshold $H_{0}$ under which outflow is zero [m] -+ `storfunc` - Type of lake storage curve ; 1 for $S = AH$ (default) and 2 for $S = - f(H)$ from lake data and interpolation ++ `storfunc` - Type of lake storage curve ; 1 for $S = AH$ (default) and 2 for $S = f(H)$ from + lake data and interpolation + `outflowfunc` - Type of lake rating curve ; 1 for $Q = f(H)$ from lake data and interpolation, 2 for general $Q = b(H - H_{0})^{e}$ and 3 in the case of Puls Approach $Q = b(H - H_{0})^{2}$ (default) + `b` - Rating curve coefficient + `e` - Rating curve exponent -By default, the lakes are not included in the model. To include them, put the following line -in the TOML file of the model: +By default, the lakes are not included in the model. To include them, put the following line in +the TOML file of the model: ```toml [model] @@ -155,8 +158,8 @@ supplied in the same folder of the TOML file. Naming of the files uses the ID of where data are available and is of the form `lake_sh_1.csv` and `lake_hq_1.csv` for respectively the storage and rating curves of lake with ID 1. -The storage curve is stored in a CSV file with lake level $\SIb{}{m}$ in the first column `H` and -corresponding lake storage $\SIb{}{m^3}$ in the second column `S`: +The storage curve is stored in a CSV file with lake level $\SIb{}{m}$ in the first column `H` +and corresponding lake storage $\SIb{}{m^3}$ in the second column `S`: ``` H, S @@ -166,12 +169,12 @@ H, S 394.21, 869719000 ``` -The rating curve uses level and discharge data depending on the Julian day of the year -(JDOY), and can be also used for regulated lakes/ dams. The first line contains `H` for the -first column. The other lines contain the water level and the corresponding discharges for -the different JDOY (1-365), see also the example below, that shows part of a CSV file (first -4 Julian days). The volume above the maximum water level of the rating curve is assumed to -flow instantaneously out of the lake (overflow). +The rating curve uses level and discharge data depending on the Julian day of the year (JDOY), +and can be also used for regulated lakes/ dams. The first line contains `H` for the first +column. The other lines contain the water level and the corresponding discharges for the +different JDOY (1-365), see also the example below, that shows part of a CSV file (first 4 +Julian days). The volume above the maximum water level of the rating curve is assumed to flow +instantaneously out of the lake (overflow). ``` H @@ -183,18 +186,17 @@ H 394.05, 52.179, 52.179, 52.179, 52.179 ``` Linked lakes: In some cases, lakes can be linked and return flow can be allowed from the -downstream to the upstream lake. The linked lakes are defined in the `linkedlakelocs` -parameter that represent the downstream lake location ID, at the grid cell of the upstream -lake location. +downstream to the upstream lake. The linked lakes are defined in the `linkedlakelocs` parameter +that represent the downstream lake location ID, at the grid cell of the upstream lake location. ::: {.callout-note} - In every file, level units are meters [m] above lake bottom and not meters above sea - level [m asl]. Especially with storage/rating curves coming from data, please be careful - and convert units if needed. +In every file, level units are meters [m] above lake bottom and not meters above sea level [m +asl]. Especially with storage/rating curves coming from data, please be careful and convert +units if needed. ::: ## References + Bos M.G., 1989. Discharge measurement structures. Third revised edition, International Institute for Land Reclamation and Improvement ILRI, Wageningen, The Netherlands. -+ Burek P., Van der Knijf J.M., Ad de Roo, 2013. LISFLOOD – Distributed Water Balance and - flood Simulation Model – Revised User Manual. DOI: http://dx.doi.org/10.2788/24719. \ No newline at end of file ++ Burek P., Van der Knijf J.M., Ad de Roo, 2013. LISFLOOD – Distributed Water Balance and flood + Simulation Model – Revised User Manual. DOI: http://dx.doi.org/10.2788/24719. \ No newline at end of file diff --git a/docs/model_docs/model_configurations.qmd b/docs/model_docs/model_configurations.qmd index 228e93dfe..148194413 100644 --- a/docs/model_docs/model_configurations.qmd +++ b/docs/model_docs/model_configurations.qmd @@ -2,12 +2,12 @@ title: Model configurations --- -There are several model configurations supported by wflow. These model configurations -require slightly different input requirements, yet the general structure is similar for each -model. A wflow model configuration consists of a `vertical` [SBM](@ref vert_sbm) concept in -combination with `lateral` concepts that control how water is routed for example over the -land or river domain. For the wflow\_sbm model different model configurations are possible. -The following model configurations are supported in wflow: +There are several model configurations supported by wflow. These model configurations require +slightly different input requirements, yet the general structure is similar for each model. A +wflow model configuration consists of a `vertical` [SBM](@ref vert_sbm) concept in combination +with `lateral` concepts that control how water is routed for example over the land or river +domain. For the wflow\_sbm model different model configurations are possible. The following +model configurations are supported in wflow: - wflow\_sbm: - SBM + kinematic wave for subsurface and surface flow @@ -17,19 +17,18 @@ The following model configurations are supported in wflow: - SBM + groundwater flow + kinematic wave for surface flow - wflow\_sediment as post processing of wflow\_sbm output -Below, some explanation will be given on how to prepare a basic wflow\_sbm -model. Example data for other model configurations is provided in the section with [sample -data](@ref sample_data). +Below, some explanation will be given on how to prepare a basic wflow\_sbm model. Example data +for other model configurations is provided in the section with [sample data](@ref sample_data). ## wflow\_sbm Wflow\_sbm represents hydrological models derived from the CQflow model (Köhler et al., 2006) -that have the [SBM](@ref vert_sbm) vertical concept in common, but can have different -lateral concepts that control how water is routed for example over the land or river domain. -The soil part of SBM is largely based on the Topog\_SBM model but has had considerable -changes over time. Topog\_SBM is specifically designed to simulate fast runoff processes in -small catchments while the wflow\_sbm model can be applied more widely. The main differences are -for the vertical concept SBM of wflow\_sbm: +that have the [SBM](@ref vert_sbm) vertical concept in common, but can have different lateral +concepts that control how water is routed for example over the land or river domain. The soil +part of SBM is largely based on the Topog\_SBM model but has had considerable changes over +time. Topog\_SBM is specifically designed to simulate fast runoff processes in small catchments +while the wflow\_sbm model can be applied more widely. The main differences are for the +vertical concept SBM of wflow\_sbm: - The unsaturated zone can be split-up in different layers - The addition of evapotranspiration losses @@ -47,13 +46,13 @@ configurations: The vertical SBM concept is explained in more detail in the following section [SBM vertical concept](@ref vert_sbm). -Topog\_SBM uses an element network based on contour lines and trajectories for water -routing. Wflow\_sbm models differ in how the lateral components river, land, and subsurface -are solved. Below the different wflow\_sbm model configurations are described. +Topog\_SBM uses an element network based on contour lines and trajectories for water routing. +Wflow\_sbm models differ in how the lateral components river, land, and subsurface are solved. +Below the different wflow\_sbm model configurations are described. ### SBM + Kinematic wave -For the lateral components of this wflow\_sbm model water is routed over a D8 network, and -the kinematic wave approach is used for river, overland and lateral subsurface flow. This is +For the lateral components of this wflow\_sbm model water is routed over a D8 network, and the +kinematic wave approach is used for river, overland and lateral subsurface flow. This is described in more detail in the section [Kinematic wave](@ref kin_wave). An overview of the different processes and fluxes in the wflow_sbm model with the kinematic @@ -61,10 +60,10 @@ wave approach for river, overland and lateral subsurface flow: ![Conceptual overview of the wflow_sbm model](../images/wflow_sbm_soil.png) -Below the mapping for this wflow\_sbm model (type `sbm`) to the vertical SBM concept -(instance of `struct SBM`) and the different lateral concepts is presented. For an -explanation about the type parameters between curly braces after the `struct` name see the -section on the model parameters. +Below the mapping for this wflow\_sbm model (type `sbm`) to the vertical SBM concept (instance +of `struct SBM`) and the different lateral concepts is presented. For an explanation about the +type parameters between curly braces after the `struct` name see the section on the model +parameters. ```julia vertical => struct SBM{T,N,M} @@ -97,7 +96,8 @@ gwf_f.value = 3.0 ``` Below the mapping for this wflow\_sbm model (type `sbm_gwf`) to the vertical SBM concept (instance of `struct SBM`) and the different lateral concepts. For an explanation about the -type parameters between curly braces after the `struct` name see the section on model parameters. +type parameters between curly braces after the `struct` name see the section on model +parameters. ```julia vertical => struct SBM{T,N,M} @@ -112,10 +112,10 @@ lateral.river.reservoir => struct SimpleReservoir{T} # optional ``` ### SBM + Local inertial river -By default the model types `sbm` and `sbm_gwf` use the kinematic wave approach for river -flow. There is also the option to use the local inertial model for river flow with an -optional 1D floodplain schematization (routing is done separately for the river channel and -floodplain), by providing the following in the TOML file: +By default the model types `sbm` and `sbm_gwf` use the kinematic wave approach for river flow. +There is also the option to use the local inertial model for river flow with an optional 1D +floodplain schematization (routing is done separately for the river channel and floodplain), by +providing the following in the TOML file: ```toml [model] @@ -123,8 +123,8 @@ river_routing = "local-inertial" # optional, default is "kinematic-wave" floodplain_1d = true # optional, default is false ``` -Only the mapping for the river component changes, as shown below. For an explanation about -the type parameters between curly braces after the `struct` name see the section on the model +Only the mapping for the river component changes, as shown below. For an explanation about the +type parameters between curly braces after the `struct` name see the section on the model parameters. ```julia @@ -141,9 +141,9 @@ overland flow, by providing the following in the TOML file: river_routing = "local-inertial" land_routing = "local-inertial" ``` -The mapping for the river and land component changes, as shown below. For an explanation -about the type parameters between curly braces after the `struct` name see the section on -the model parameters. +The mapping for the river and land component changes, as shown below. For an explanation about +the type parameters between curly braces after the `struct` name see the section on the model +parameters. ```julia lateral.river => struct ShallowWaterRiver{T,R,L} @@ -159,31 +159,31 @@ catchment level are intricately linked to the processes governing sediment dynam nutrients such as phosphorus, carbon or other pollutants such as metals are influenced by sediment properties in processes such as mobilization, flocculation or deposition. To better assert and model water quality in inland systems, a better comprehension and modelling of -sediment sources and fate in the river is needed at a spatial and time scale relevant to -such issues. +sediment sources and fate in the river is needed at a spatial and time scale relevant to such +issues. The wflow\_sediment model was developed to answer such issues. It is a distributed physics-based model, based on the distributed hydrologic wflow\_sbm model. It is able to -simulate both land and in-stream processes, and relies on available global datasets, -parameter estimation and small calibration effort. +simulate both land and in-stream processes, and relies on available global datasets, parameter +estimation and small calibration effort. In order to model the exports of terrestrial sediment to the coast through the Land Ocean Aquatic Continuum or LOAC (inland waters network such as streams, lakes...), two different modelling parts were considered. The first part, called the inland sediment model, is the -modelling and estimation of soil loss and sediment yield to the river system by land -erosion, separated into vertical [Soil Erosion](@ref) processes and lateral [Sediment Flux -in overland flow](@ref). The second part, called the [River Sediment Model](@ref) is the -transport and processes of the sediment in the river system. The two parts together -constitute the wflow\_sediment model. +modelling and estimation of soil loss and sediment yield to the river system by land erosion, +separated into vertical [Soil Erosion](@ref) processes and lateral [Sediment Flux in overland +flow](@ref). The second part, called the [River Sediment Model](@ref) is the transport and +processes of the sediment in the river system. The two parts together constitute the +wflow\_sediment model. Overview of the concepts of the wflow\_sediment model: ![wflow_sediment](../images/wflow_sediment.png) ### Configuration -As sediment generation and transport processes are linked to the hydrology and water flows, -the inputs to the wflow\_sediment model come directly from a hydrological model. The -required dynamic inputs to run wflow\_sediment are: +As sediment generation and transport processes are linked to the hydrology and water flows, the +inputs to the wflow\_sediment model come directly from a hydrological model. The required +dynamic inputs to run wflow\_sediment are: - Precipitation (can also come from the hydrological forcing data), - Land runoff (overland flow) from the kinematic wave, diff --git a/docs/model_docs/parameters_lateral.qmd b/docs/model_docs/parameters_lateral.qmd index 94c97df27..f97607965 100644 --- a/docs/model_docs/parameters_lateral.qmd +++ b/docs/model_docs/parameters_lateral.qmd @@ -5,13 +5,13 @@ title: Lateral concepts ## Kinematic wave ### Surface flow -The Table below shows the parameters (fields) of struct `SurfaceFlowRiver` used for river -flow, including a description of these parameters, the unit, and default value if -applicable. The parameters in bold represent model parameters that can be set through static -input data (netCDF), and can be listed in the TOML configuration file under -`[input.lateral.river]` to map the internal model parameter to the external netCDF variable. -The input parameter `slope` (listed under `[input.lateral.river]`) is not equal to the -internal model parameter `sl`, and is listed in the Table below between parentheses. +The Table below shows the parameters (fields) of struct `SurfaceFlowRiver` used for river flow, +including a description of these parameters, the unit, and default value if applicable. The +parameters in bold represent model parameters that can be set through static input data +(netCDF), and can be listed in the TOML configuration file under `[input.lateral.river]` to map +the internal model parameter to the external netCDF variable. The input parameter `slope` +(listed under `[input.lateral.river]`) is not equal to the internal model parameter `sl`, and +is listed in the Table below between parentheses. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -47,12 +47,12 @@ internal model parameter `sl`, and is listed in the Table below between parenthe : {.striped .hover} The Table below shows the parameters (fields) of struct `SurfaceFlowLand` used for overland -flow, including a description of these parameters, the unit, and default value if -applicable. The parameters in bold represent model parameters that can be set through static -input data (netCDF), and can be listed in the TOML configuration file under -`[input.lateral.land]` to map the internal model parameter to the external netCDF variable. -The input parameter `slope` (listed under `[input.lateral.land]`) is not equal to the -internal model parameter `sl`, and is listed in the Table below between parentheses. +flow, including a description of these parameters, the unit, and default value if applicable. +The parameters in bold represent model parameters that can be set through static input data +(netCDF), and can be listed in the TOML configuration file under `[input.lateral.land]` to map +the internal model parameter to the external netCDF variable. The input parameter `slope` +(listed under `[input.lateral.land]`) is not equal to the internal model parameter `sl`, and is +listed in the Table below between parentheses. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -80,41 +80,41 @@ internal model parameter `sl`, and is listed in the Table below between parenthe : {.striped .hover} ### Lateral subsurface flow -The Table below shows the parameters (fields) of struct `LateralSSF`, including a -description of these parameters, the unit, and default value if applicable. The parameters -in bold represent model parameters that can be set through static input data (netCDF). The -soil related parameters `f`, `soilthickness`, `z_exp`, `theta_s` and `theta_r` are derived from the -vertical `SBM` concept (including unit conversion for `f`, `z_exp` and `soilthickness`), and -can be listed in the TOML configuration file under `[input.vertical]`, to map the internal -model parameter to the external netCDF variable. The internal slope model parameter `slope` is -set through the TOML file as follows: +The Table below shows the parameters (fields) of struct `LateralSSF`, including a description +of these parameters, the unit, and default value if applicable. The parameters in bold +represent model parameters that can be set through static input data (netCDF). The soil related +parameters `f`, `soilthickness`, `z_exp`, `theta_s` and `theta_r` are derived from the vertical +`SBM` concept (including unit conversion for `f`, `z_exp` and `soilthickness`), and can be +listed in the TOML configuration file under `[input.vertical]`, to map the internal model +parameter to the external netCDF variable. The internal slope model parameter `slope` is set +through the TOML file as follows: ```toml [input.lateral.land] slope = "Slope" ``` -The parameter `kh_0` is computed by multiplying the vertical hydraulic conductivity at the -soil surface `kv_0` (including unit conversion) of the vertical `SBM` concept with the -internal parameter `khfrac` \[-\] (default value of 1.0). The internal model parameter -`khfrac` is set through the TOML file as follows: +The parameter `kh_0` is computed by multiplying the vertical hydraulic conductivity at the soil +surface `kv_0` (including unit conversion) of the vertical `SBM` concept with the internal +parameter `khfrac` \[-\] (default value of 1.0). The internal model parameter `khfrac` is set +through the TOML file as follows: ```toml [input.lateral.subsurface] ksathorfrac = "KsatHorFrac" ``` -The `khfrac` parameter compensates for anisotropy, small scale `kv_0` measurements (soil -core) that do not represent larger scale hydraulic conductivity, and smaller flow length -scales (hillslope) in reality, not represented by the model resolution. +The `khfrac` parameter compensates for anisotropy, small scale `kv_0` measurements (soil core) +that do not represent larger scale hydraulic conductivity, and smaller flow length scales +(hillslope) in reality, not represented by the model resolution. -For the vertical [SBM](@ref params_sbm) concept different vertical hydraulic conductivity -depth profiles are possible, and these also determine which `LateralSSF` parameters are used +For the vertical [SBM](@ref params_sbm) concept different vertical hydraulic conductivity depth +profiles are possible, and these also determine which `LateralSSF` parameters are used including the input requirements for the computation of lateral subsurface flow. For the `exponential` profile the model parameters `kh_0` and `f` are used. For the `exponential_constant` profile `kh_0` and `f` are used, and `z_exp` is required as part of -`[input.vertical]`. For the `layered` profile, `SBM` model parameter `kv` is used, and for -the `layered_exponential` profile `kv` is used and `z_exp` is required as part of +`[input.vertical]`. For the `layered` profile, `SBM` model parameter `kv` is used, and for the +`layered_exponential` profile `kv` is used and `z_exp` is required as part of `[input.vertical]`. | Parameter | Description | Unit | Default | @@ -144,11 +144,11 @@ the `layered_exponential` profile `kv` is used and `z_exp` is required as part o ### River flow The Table below shows the parameters (fields) of struct `ShallowWaterRiver`, including a -description of these parameters, the unit, and default value if applicable. The parameters -in bold represent model parameters that can be set through static input data (netCDF), and -can be listed in the TOML configuration file under `[input.lateral.river]`, to map the -internal model parameter to the external netCDF variable. The parameter river bed elevation -`zb` is based on the bankfull elevation and depth input data: +description of these parameters, the unit, and default value if applicable. The parameters in +bold represent model parameters that can be set through static input data (netCDF), and can be +listed in the TOML configuration file under `[input.lateral.river]`, to map the internal model +parameter to the external netCDF variable. The parameter river bed elevation `zb` is based on +the bankfull elevation and depth input data: ```toml [input.lateral.river] @@ -159,8 +159,8 @@ bankfull_depth = "RiverDepth" When floodplain routing (parameter `floodplain`) is included as part of local inertial river flow, parameter `q_av` represents the total average discharge of the river channel and floodplain routing, and parameter `q_channel_av` represents average river channel discharge. -Otherwise parameters `q_av` and `q_channel_av` represent both average river channel -discharge (are equal). +Otherwise parameters `q_av` and `q_channel_av` represent both average river channel discharge +(are equal). The input parameter `n` (listed under `[input.lateral.river]`) is not equal to the internal model parameter `mannings_n`, and is listed in the Table below between parentheses. @@ -216,13 +216,13 @@ model parameter `mannings_n`, and is listed in the Table below between parenthes ### 1D floodplain The Table below shows the parameters (fields) of struct `FloodPlain` (part of struct -`ShallowWaterRiver`), including a description of these parameters, the unit, and default -value if applicable. The parameters in bold represent model parameters that can be set -through static input data (netCDF), and can be listed in the TOML configuration file under -`[input.lateral.river.floodplain]`, to map the internal model parameter to the external -netCDF variable. The input parameter `n` (listed under `[input.lateral.river.floodplain]`) -is not equal to the internal model parameter `mannings_n`, and is listed in the Table below -between parentheses. +`ShallowWaterRiver`), including a description of these parameters, the unit, and default value +if applicable. The parameters in bold represent model parameters that can be set through static +input data (netCDF), and can be listed in the TOML configuration file under +`[input.lateral.river.floodplain]`, to map the internal model parameter to the external netCDF +variable. The input parameter `n` (listed under `[input.lateral.river.floodplain]`) is not +equal to the internal model parameter `mannings_n`, and is listed in the Table below between +parentheses. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -267,13 +267,13 @@ internal model parameter `depth`, and is listed in the Table below between paren ### Overland flow The Table below shows the parameters (fields) of struct `ShallowWaterLand`, including a -description of these parameters, the unit, and default value if applicable. The parameters -in bold represent model parameters that can be set through static input data (netCDF), and -can be listed in the TOML configuration file under `[input.lateral.land]`, to map the -internal model parameter to the external netCDF variable. +description of these parameters, the unit, and default value if applicable. The parameters in +bold represent model parameters that can be set through static input data (netCDF), and can be +listed in the TOML configuration file under `[input.lateral.land]`, to map the internal model +parameter to the external netCDF variable. -The mannings roughness (for the computation of `mannings_n_sq`) should be provided as -follows in the TOML file: +The mannings roughness (for the computation of `mannings_n_sq`) should be provided as follows +in the TOML file: ```toml [input.lateral.land] @@ -313,8 +313,8 @@ internal model parameter `z`, and is listed in the Table below between parenthes ## Water allocation river The Table below shows the parameters (fields) of struct `AllocationRiver`, used when water -demand and allocation is computed (optional), including a description of these parameters, -the unit, and default value if applicable. +demand and allocation is computed (optional), including a description of these parameters, the +unit, and default value if applicable. | parameter | description | unit | default | |:--------------- | ------------------| ----- | -------- | @@ -344,14 +344,14 @@ description of these parameters, the unit, and default value if applicable. Stru ### Unconfined aquifer The Table below shows the parameters (fields) of struct `UnconfinedAquifer`, including a -description of these parameters, the unit, and default value if applicable. The parameters -in bold represent model parameters that can be set through static input data (netCDF), and -can be listed in the TOML configuration file under `[input.lateral.subsurface]`, to map the -internal model parameter to the external netCDF variable. For some input parameters the -parameter listed under `[input.lateral.subsurface]` is not equal to the internal model -parameter, these are listed in the Table below between parentheses after the internal model -parameter. The `top` parameter is provided by the external parameter `altitude` as part of -the static input data and set as follows through the TOML file: +description of these parameters, the unit, and default value if applicable. The parameters in +bold represent model parameters that can be set through static input data (netCDF), and can be +listed in the TOML configuration file under `[input.lateral.subsurface]`, to map the internal +model parameter to the external netCDF variable. For some input parameters the parameter listed +under `[input.lateral.subsurface]` is not equal to the internal model parameter, these are +listed in the Table below between parentheses after the internal model parameter. The `top` +parameter is provided by the external parameter `altitude` as part of the static input data and +set as follows through the TOML file: ```toml [input] @@ -359,8 +359,8 @@ the static input data and set as follows through the TOML file: altitude = "wflow_dem" ``` -The input parameter `conductivity` (listed under `[input.lateral.subsurface]`) is not equal -to the internal model parameter `kh_0`, and is listed in the Table below between parentheses. +The input parameter `conductivity` (listed under `[input.lateral.subsurface]`) is not equal to +the internal model parameter `kh_0`, and is listed in the Table below between parentheses. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -375,13 +375,13 @@ to the internal model parameter `kh_0`, and is listed in the Table below between : {.striped .hover} ### Constant Head -The Table below shows the parameters (fields) of struct `ConstantHead`, including a -description of these parameters, the unit, and default value if applicable. The parameters -in bold represent model parameters that can be set through static input data (netCDF), and -can be listed in the TOML configuration file under `[input.lateral.subsurface]`, to map the -internal model parameter to the external netCDF variable. The input parameter -`constant_head` (listed under `[input.lateral.subsurface]`) is not equal to the internal -model parameter `head`, and is listed in the Table below between parentheses. +The Table below shows the parameters (fields) of struct `ConstantHead`, including a description +of these parameters, the unit, and default value if applicable. The parameters in bold +represent model parameters that can be set through static input data (netCDF), and can be +listed in the TOML configuration file under `[input.lateral.subsurface]`, to map the internal +model parameter to the external netCDF variable. The input parameter `constant_head` (listed +under `[input.lateral.subsurface]`) is not equal to the internal model parameter `head`, and is +listed in the Table below between parentheses. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -393,12 +393,12 @@ model parameter `head`, and is listed in the Table below between parentheses. #### River The Table below shows the parameters (fields) of struct `River`, including a description of -these parameters, the unit, and default value if applicable. The parameters in bold -represent model parameters that can be set through static input data (netCDF), and can be -listed in the TOML configuration file under `[input.lateral.subsurface]`, to map the -internal model parameter to the external netCDF variable. The input parameter `river_bottom` -(listed under `[input.lateral.subsurface]`) is not equal to the internal model parameter -`bottom`, and is listed in the Table below between parentheses. +these parameters, the unit, and default value if applicable. The parameters in bold represent +model parameters that can be set through static input data (netCDF), and can be listed in the +TOML configuration file under `[input.lateral.subsurface]`, to map the internal model parameter +to the external netCDF variable. The input parameter `river_bottom` (listed under +`[input.lateral.subsurface]`) is not equal to the internal model parameter `bottom`, and is +listed in the Table below between parentheses. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -411,14 +411,13 @@ internal model parameter to the external netCDF variable. The input parameter `r : {.striped .hover} #### Drainage -The Table below shows the parameters (fields) of struct `Drainage`, including a description -of these parameters, the unit, and default value if applicable. The parameters in bold -represent model parameters that can be set through static input data (netCDF), and can be -listed in the TOML configuration file under `[input.lateral.subsurface]`, to map the -internal model parameter to the external netCDF variable. For some input parameters the -parameter listed under `[input.lateral.subsurface]` is not equal to the internal model -parameter, these are listed in the Table below between parentheses after the internal model -parameter. +The Table below shows the parameters (fields) of struct `Drainage`, including a description of +these parameters, the unit, and default value if applicable. The parameters in bold represent +model parameters that can be set through static input data (netCDF), and can be listed in the +TOML configuration file under `[input.lateral.subsurface]`, to map the internal model parameter +to the external netCDF variable. For some input parameters the parameter listed under +`[input.lateral.subsurface]` is not equal to the internal model parameter, these are listed in +the Table below between parentheses after the internal model parameter. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -429,8 +428,8 @@ parameter. : {.striped .hover} #### Recharge -The Table below shows the parameters (fields) of struct `Recharge`, including a description -of these parameters, the unit, and default value if applicable. +The Table below shows the parameters (fields) of struct `Recharge`, including a description of +these parameters, the unit, and default value if applicable. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -440,8 +439,8 @@ of these parameters, the unit, and default value if applicable. : {.striped .hover} #### Head boundary -The Table below shows the parameters (fields) of struct `HeadBoundary`, including a -description of these parameters, the unit, and default value if applicable. +The Table below shows the parameters (fields) of struct `HeadBoundary`, including a description +of these parameters, the unit, and default value if applicable. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -493,10 +492,10 @@ description of these parameters, the unit, and default value if applicable. : {.striped .hover} ### River flow -The Table below shows external parameters that can be set through static input data -(netCDF), and can be listed in the TOML configuration file under `[input.lateral.river]`. -These external parameters are not part of struct `RiverSediment`, but used to calculate -parameters of struct `RiverSediment`. +The Table below shows external parameters that can be set through static input data (netCDF), +and can be listed in the TOML configuration file under `[input.lateral.river]`. These external +parameters are not part of struct `RiverSediment`, but used to calculate parameters of struct +`RiverSediment`. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -510,14 +509,14 @@ parameters of struct `RiverSediment`. : {.striped .hover} The Table below shows the parameters (fields) of struct `RiverSediment`, including a -description of these parameters, the unit, and default value if applicable. The parameters -in bold represent model parameters that can be set through static and forcing input data -(netCDF), and can be listed in the TOML configuration file under `[input.lateral.river]`, to -map the internal model parameter to the external netCDF variable. For some input parameters -the parameter listed under `[input.lateral.river]` is not equal to the internal model -parameter, these are listed in the Table below between parentheses after the internal model -parameter. For example, internal model parameter `sl` is mapped as follows in the TOML file -to the external netCDF variable `RiverSlope`: +description of these parameters, the unit, and default value if applicable. The parameters in +bold represent model parameters that can be set through static and forcing input data (netCDF), +and can be listed in the TOML configuration file under `[input.lateral.river]`, to map the +internal model parameter to the external netCDF variable. For some input parameters the +parameter listed under `[input.lateral.river]` is not equal to the internal model parameter, +these are listed in the Table below between parentheses after the internal model parameter. For +example, internal model parameter `sl` is mapped as follows in the TOML file to the external +netCDF variable `RiverSlope`: ```toml [input.vertical] diff --git a/docs/model_docs/parameters_vertical.qmd b/docs/model_docs/parameters_vertical.qmd index 792bef7de..4823dae03 100644 --- a/docs/model_docs/parameters_vertical.qmd +++ b/docs/model_docs/parameters_vertical.qmd @@ -3,24 +3,23 @@ title: Vertical concepts --- ## SBM -The Table below shows the parameters (fields) of struct `SBM`, including a description of -these parameters, the unit, and default value if applicable. The parameters in bold -represent model parameters that can be set through static and forcing input data (netCDF), -and can be listed in the TOML configuration file under `[input.vertical]`, to map the -internal model parameter to the external netCDF variable. For some input parameters the -parameter listed under `[input.vertical]` is not equal to the internal model parameter, -these are listed in the Table below between parentheses after the internal model parameter. -For example, internal model parameter `sl` is mapped as follows in the TOML file to the -external netCDF variable `Sl`: +The Table below shows the parameters (fields) of struct `SBM`, including a description of these +parameters, the unit, and default value if applicable. The parameters in bold represent model +parameters that can be set through static and forcing input data (netCDF), and can be listed in +the TOML configuration file under `[input.vertical]`, to map the internal model parameter to +the external netCDF variable. For some input parameters the parameter listed under +`[input.vertical]` is not equal to the internal model parameter, these are listed in the Table +below between parentheses after the internal model parameter. For example, internal model +parameter `sl` is mapped as follows in the TOML file to the external netCDF variable `Sl`: ```toml [input.vertical] specific_leaf = "Sl" ``` -Different [vertical hydraulic conductivity depth profiles](@ref soil): `exponential` -(default), `exponential_constant`, `layered` and `layered_exponential` can be provided -through the TOML file. Below an example for the `exponential_constant` profile: +Different [vertical hydraulic conductivity depth profiles](@ref soil): `exponential` (default), +`exponential_constant`, `layered` and `layered_exponential` can be provided through the TOML +file. Below an example for the `exponential_constant` profile: ```toml [input.vertical] @@ -29,8 +28,8 @@ ksat_profile = "exponential_constant" For the `exponential` profile the input parameters `kv_0` and `f` are used. For the `exponential_constant` profile `kv_0` and `f` are used, and `z_exp` is required as input. For -the `layered` profile, input parameter `kv` is used, and for the `layered_exponential` -profile `kv` is used and `z_layered` is required as input. +the `layered` profile, input parameter `kv` is used, and for the `layered_exponential` profile +`kv` is used and `z_layered` is required as input. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -155,10 +154,10 @@ profile `kv` is used and `z_layered` is required as input. ### Paddy The Table below shows the parameters (fields) of struct `Paddy`, including a description of -these parameters, the unit, and default value if applicable. The parameters in bold -represent model parameters that can be set through static and forcing input data (netCDF), -and can be listed in the TOML configuration file under `[input.vertical.paddy]`, to map the -internal model parameter to the external netCDF variable. +these parameters, the unit, and default value if applicable. The parameters in bold represent +model parameters that can be set through static and forcing input data (netCDF), and can be +listed in the TOML configuration file under `[input.vertical.paddy]`, to map the internal model +parameter to the external netCDF variable. | parameter | description | unit | default | |:---------------| --------------- | ---------------------- | ----- | @@ -173,11 +172,11 @@ internal model parameter to the external netCDF variable. | `h` | actual water depth in paddy field | mm | - | ### Non-paddy -The Table below shows the parameters (fields) of struct `NonPaddy`, including a description -of these parameters, the unit, and default value if applicable. The parameters in bold -represent model parameters that can be set through static and forcing input data (netCDF), -and can be listed in the TOML configuration file under `[input.vertical.nonpaddy]`, to map -the internal model parameter to the external netCDF variable. +The Table below shows the parameters (fields) of struct `NonPaddy`, including a description of +these parameters, the unit, and default value if applicable. The parameters in bold represent +model parameters that can be set through static and forcing input data (netCDF), and can be +listed in the TOML configuration file under `[input.vertical.nonpaddy]`, to map the internal +model parameter to the external netCDF variable. | parameter | description | unit | default | |:---------------| --------------- | ---------------------- | ----- | @@ -189,12 +188,12 @@ the internal model parameter to the external netCDF variable. ### Non-irrigation (industry, domestic and livestock) The Table below shows the parameters (fields) of struct `NonIrrigationDemand`, including a -description of these parameters, the unit, and default value if applicable. The parameters -in bold represent model parameters that can be set through static and forcing input data -(netCDF). These parameters can be listed for the sectors industry, domestic and livestock, -in the TOML configuration file under `[input.vertical.industry]`, -`[input.vertical.domestic]` and `[input.vertical.livestock]`, to map the internal model -parameter to the external netCDF variable. +description of these parameters, the unit, and default value if applicable. The parameters in +bold represent model parameters that can be set through static and forcing input data (netCDF). +These parameters can be listed for the sectors industry, domestic and livestock, in the TOML +configuration file under `[input.vertical.industry]`, `[input.vertical.domestic]` and +`[input.vertical.livestock]`, to map the internal model parameter to the external netCDF +variable. | parameter | description | unit | default | |:---------------| --------------- | ---------------------- | ----- | @@ -205,11 +204,10 @@ parameter to the external netCDF variable. ### Water allocation land The Table below shows the parameters (fields) of struct `AllocationLand`, including a -description of these parameters, the unit, and default value if applicable. The parameters -in bold represent model parameters that can be set through static and forcing input data -(netCDF), and can be listed in the TOML configuration file under -`[input.vertical.allocation]`, to map the internal model parameter to the external netCDF -variable. +description of these parameters, the unit, and default value if applicable. The parameters in +bold represent model parameters that can be set through static and forcing input data (netCDF), +and can be listed in the TOML configuration file under `[input.vertical.allocation]`, to map +the internal model parameter to the external netCDF variable. | parameter | description | unit | default | |:---------------| --------------- | ---------------------- | ----- | @@ -232,10 +230,10 @@ variable. ## Sediment -The Table below shows external parameters that can be set through static input data -(netCDF), and can be listed in the TOML configuration file under `[input.vertical]`. These -external parameters are not part of struct `LandSediment`, but used to calculate parameters -of struct `LandSediment`. +The Table below shows external parameters that can be set through static input data (netCDF), +and can be listed in the TOML configuration file under `[input.vertical]`. These external +parameters are not part of struct `LandSediment`, but used to calculate parameters of struct +`LandSediment`. | Parameter | Description | Unit | Default | |:----| -------- | ---- | ---- | @@ -245,15 +243,14 @@ of struct `LandSediment`. | `lakeareas` | lake coverage | - | - | : {.striped .hover} -The Table below shows the parameters (fields) of struct `LandSediment`, including a -description of these parameters, the unit, and default value if applicable. The parameters -in bold represent model parameters that can be set through static and forcing input data -(netCDF), and can be listed in the TOML configuration file under `[input.vertical]`, to map -the internal model parameter to the external netCDF variable. For some input parameters the -parameter listed under `[input.vertical]` is not equal to the internal model parameter, -these are listed in the Table below between parentheses after the internal model parameter. -For example, internal model parameter `sl` is mapped as follows in the TOML file to the -external netCDF variable `Sl`: +The Table below shows the parameters (fields) of struct `LandSediment`, including a description +of these parameters, the unit, and default value if applicable. The parameters in bold +represent model parameters that can be set through static and forcing input data (netCDF), and +can be listed in the TOML configuration file under `[input.vertical]`, to map the internal +model parameter to the external netCDF variable. For some input parameters the parameter listed +under `[input.vertical]` is not equal to the internal model parameter, these are listed in the +Table below between parentheses after the internal model parameter. For example, internal model +parameter `sl` is mapped as follows in the TOML file to the external netCDF variable `Sl`: ```toml [input.vertical] diff --git a/docs/model_docs/vertical/sbm.qmd b/docs/model_docs/vertical/sbm.qmd index 6c9feab0e..e553c9bd1 100644 --- a/docs/model_docs/vertical/sbm.qmd +++ b/docs/model_docs/vertical/sbm.qmd @@ -19,9 +19,9 @@ performed based on the air temperature. If the temperature is below a threshold (`tt`), precipitation will fall as snow. An interval parameter (`tti`) defines the range over which precipitation is partly falling as snow, and partly as rain. Snowfall is added to the snowpack, where it is subject to melting and refreezing (see the section on [snow and -glaciers](@ref snow)). The amount of rainfall is subject to [interception](@ref -interception), and ultimately becomes available for [evaporation](@ref evap) and/or [soil -processes](@ref soil). +glaciers](@ref snow)). The amount of rainfall is subject to [interception](@ref interception), +and ultimately becomes available for [evaporation](@ref evap) and/or [soil processes](@ref +soil). ![Division between snow and precipitation based on the threshold temperature](../../images/snowfall.png) @@ -35,15 +35,15 @@ processes](@ref soil). ## Rainfall interception Two different interception models are available: the analytical Gash model, and the modified -Rutter model. The simulation timestep defines which interception model is used, where daily -(or larger) timesteps use the Gash model, and timesteps smaller than daily use the modified -Rutter model. +Rutter model. The simulation timestep defines which interception model is used, where daily (or +larger) timesteps use the Gash model, and timesteps smaller than daily use the modified Rutter +model. ### The analytical (Gash) model (Gash, 1979) -The analytical model of rainfall interception is based on Rutter's numerical model. Simplifications -allow the model to be applied on a daily basis, although a -storm-based approach will yield better results in situations with more than one storm per -day. The amount of water needed to completely saturate the canopy is defined as: +The analytical model of rainfall interception is based on Rutter's numerical model. +Simplifications allow the model to be applied on a daily basis, although a storm-based approach +will yield better results in situations with more than one storm per day. The amount of water +needed to completely saturate the canopy is defined as: $$ P'=\frac{-\overline{R}S}{\overline{E}_{w}}\log\left[1-\frac{\overline{E}_{w}}{\overline{R}}(1-p-p_{t})^{-1}\right] @@ -51,14 +51,13 @@ $$ where $\overline{R}$ is the average precipitation intensity on a saturated canopy and $\overline{E}_{w}$ the average evaporation from the wet canopy and with the vegetation -parameters $S$, $p$ and $p_t$ as defined previously. The model uses a series of -expressions to calculate the interception loss during different phases of a storm. An -analytical integration of the total evaporation and rainfall under saturated canopy -conditions is performed for each storm to determine average values of $\overline{E}_{w}$ -and $\overline{R}$. The total evaporation from the canopy (the total interception loss) is -calculated as the sum of the components listed in the table below. Interception losses from -the stems are calculated for days with $P\geq S_{t}/p_{t}$. $p_t$ and $S_t$ are small -and neglected. +parameters $S$, $p$ and $p_t$ as defined previously. The model uses a series of expressions to +calculate the interception loss during different phases of a storm. An analytical integration +of the total evaporation and rainfall under saturated canopy conditions is performed for each +storm to determine average values of $\overline{E}_{w}$ and $\overline{R}$. The total +evaporation from the canopy (the total interception loss) is calculated as the sum of the +components listed in the table below. Interception losses from the stems are calculated for +days with $P\geq S_{t}/p_{t}$. $p_t$ and $S_t$ are small and neglected. Table: Formulation of the components of interception loss according to Gash: @@ -72,25 +71,25 @@ Table: Formulation of the components of interception loss according to Gash: | Evaporation from trunks in $m+n-q$ storms that do not fill the trunk storage | $p_{t}\sum_{j=1}^{m+n-q}P_{g,j}$ | In applying the analytical model, saturated conditions are assumed to occur when the hourly -rainfall exceeds a certain threshold. Often a threshold of 0.5 mm/hr is used. -$\overline{R}$ is calculated for all hours when the rainfall exceeds the threshold to give -an estimate of the mean rainfall rate onto a saturated canopy. +rainfall exceeds a certain threshold. Often a threshold of 0.5 mm/hr is used. $\overline{R}$ is +calculated for all hours when the rainfall exceeds the threshold to give an estimate of the +mean rainfall rate onto a saturated canopy. -Gash (1979) has shown that in a regression of interception loss on rainfall (on a storm -basis) the regression coefficient should equal to $\overline{E}_w/\overline{R}$. Assuming -that neither $\overline{E}_w$ nor $\overline{R}$ vary considerably in time, -$\overline{E}_w$ can be estimated in this way from $\overline{R}$ in the absence of -above-canopy climatic observations. Values derived in this way generally tend to be (much) -higher than those calculated with the penman-monteith equation. +Gash (1979) has shown that in a regression of interception loss on rainfall (on a storm basis) +the regression coefficient should equal to $\overline{E}_w/\overline{R}$. Assuming that neither +$\overline{E}_w$ nor $\overline{R}$ vary considerably in time, $\overline{E}_w$ can be +estimated in this way from $\overline{R}$ in the absence of above-canopy climatic observations. +Values derived in this way generally tend to be (much) higher than those calculated with the +penman-monteith equation. ### The modified rutter model For sub daily timesteps the interception is calculated using a simplification of the Rutter model. The simplified model is solved explicitly and does not take drainage from the canopy -into account. The amount of stemflow is taken as a fraction (`0.1 * canopygapfraction`) of -the precipitation. Throughfall equals to the amount of water that cannot be stored by the -canopy, plus the rainfall that is not captured by the canopy. Water can evaporate from the -canopy storage, taken as the minimum between potential evaporation and the current storage. -The "left-over" potential evaporation (if any) is returned as output. +into account. The amount of stemflow is taken as a fraction (`0.1 * canopygapfraction`) of the +precipitation. Throughfall equals to the amount of water that cannot be stored by the canopy, +plus the rainfall that is not captured by the canopy. Water can evaporate from the canopy +storage, taken as the minimum between potential evaporation and the current storage. The +"left-over" potential evaporation (if any) is returned as output. ```@docs Wflow.rainfall_interception_modrut @@ -98,8 +97,8 @@ Wflow.rainfall_interception_modrut ### Interception parameters from LAI The SBM concept can determine the interception parameters from leaf area index (LAI) -climatology. In order to switch this on you must define this cyclic parameter in the TOML -file, the parameter is read from `path_static`, as follows: +climatology. In order to switch this on you must define this cyclic parameter in the TOML file, +the parameter is read from `path_static`, as follows: ```toml [input] @@ -113,9 +112,9 @@ Furthermore, these additional parameters are required: + Storage woody part of vegetation (`swood` ``\SIb{}{mm}``) + Extinction coefficient (`kext` ``\SIb{}{-}``) -Here it is assumed that `cmax` ``\SIb{}{mm}`` (leaves) (canopy storage capacity for the leaves only) -relates linearly with LAI (c.f. Van Dijk and Bruijnzeel 2001). This is done via the `sl`. `sl` -can be determined through a lookup table with land cover based on literature (Pitman 1989, +Here it is assumed that `cmax` ``\SIb{}{mm}`` (leaves) (canopy storage capacity for the leaves +only) relates linearly with LAI (c.f. Van Dijk and Bruijnzeel 2001). This is done via the `sl`. +`sl` can be determined through a lookup table with land cover based on literature (Pitman 1989, Lui 1998). Next the `cmax` (leaves) is determined using: $$ @@ -145,24 +144,24 @@ wflow\_sbm is used for a surface completely covered by vegetation, and does not effect of growing stages of vegetation and soil cover. These effects are handled separately through the use of the canopy gap fraction. -It is assumed that the potential evaporation rate of intercepted water by vegetation is -equal to the potential evapotranspiration rate of vegetation (fully covering the soil) -multiplied by the canopy fraction. The potential evapotranspiration rate left over after -interception is available for transpiration. For potential open water evaporation (river and -water bodies) the potential reference evapotranspiration rate is used (multipled by the -river fraction `riverfrac`, and open water fraction `waterfrac`). Also for potential soil -evaporation the potential reference evapotranspiration rate is used, multiplied by the -canopy gap fraction corrected by the sum of total water fraction (`riverfrac` and -`waterfrac`) and the fraction covered by a glacier (`glacierfrac`). +It is assumed that the potential evaporation rate of intercepted water by vegetation is equal +to the potential evapotranspiration rate of vegetation (fully covering the soil) multiplied by +the canopy fraction. The potential evapotranspiration rate left over after interception is +available for transpiration. For potential open water evaporation (river and water bodies) the +potential reference evapotranspiration rate is used (multipled by the river fraction +`riverfrac`, and open water fraction `waterfrac`). Also for potential soil evaporation the +potential reference evapotranspiration rate is used, multiplied by the canopy gap fraction +corrected by the sum of total water fraction (`riverfrac` and `waterfrac`) and the fraction +covered by a glacier (`glacierfrac`). ### Bare soil evaporation -If there is only one soil layer present in the wflow\_sbm model, the bare soil evaporation -is scaled according to the wetness of the soil layer. The fraction of bare soil is assumed -to be equal to the fraction not covered by the canopy (`canopygapfraction`) corrected by the -total water fraction. When the soil is fully saturated, evaporation is set to equal the -potential reference evaporation. When the soil is not fully saturated, actual evaporation -decreases linearly with decreasing soil moisture values, as indicated by the figure below. +If there is only one soil layer present in the wflow\_sbm model, the bare soil evaporation is +scaled according to the wetness of the soil layer. The fraction of bare soil is assumed to be +equal to the fraction not covered by the canopy (`canopygapfraction`) corrected by the total +water fraction. When the soil is fully saturated, evaporation is set to equal the potential +reference evaporation. When the soil is not fully saturated, actual evaporation decreases +linearly with decreasing soil moisture values, as indicated by the figure below. ![Evaporation reduction as function of available soil moisture](../../images/soil_evap.png) @@ -171,38 +170,38 @@ decreases linearly with decreasing soil moisture values, as indicated by the fig # https://gist.github.com/JoostBuitink/21dd32e71fd1360117fcd1c532c4fd9d#file-sbm_soil_figs-py # hide ``` --> -When more soil layers are present, soil evaporation is only provided from the upper soil -layer, and soil evaporation is split in evaporation from the unsaturated store and -evaporation from the saturated store. Water is first evaporated from the unsaturated store. -The remaining potential soil evaporation can be used for evaporation from the saturated -store, but only when the water table is present in the upper soil layer. Both the -evaporation from the unsaturated store and the evaporation from the saturated store are -limited by the minimum of the remaining potential soil evaporation and the available water -in the unsaturated/saturated zone of the upper soil layer. Also for multiple soil layers, -the evaporation (both unsaturated and saturated) decreases linearly with decreasing water -availability. +When more soil layers are present, soil evaporation is only provided from the upper soil layer, +and soil evaporation is split in evaporation from the unsaturated store and evaporation from +the saturated store. Water is first evaporated from the unsaturated store. The remaining +potential soil evaporation can be used for evaporation from the saturated store, but only when +the water table is present in the upper soil layer. Both the evaporation from the unsaturated +store and the evaporation from the saturated store are limited by the minimum of the remaining +potential soil evaporation and the available water in the unsaturated/saturated zone of the +upper soil layer. Also for multiple soil layers, the evaporation (both unsaturated and +saturated) decreases linearly with decreasing water availability. ### Transpiration The maximum possible root water extraction rate for each soil layer is determined by -partitioning the potential transpiration rate ``T_p`` based on the fraction of the total -root length (`rootfraction` ``\SIb{}{-}``) in each soil layer. A root water uptake reduction model is -used to calculate a reduction coefficient as a function of soil water pressure, that may -reduce the maximum possible root water extraction rate. The root water uptake reduction -model is based on the concept proposed by Feddes et al. (1978). This concept defines a -reduction coefficient ``\SIb{\alpha}{-}`` as a function of soil water pressure (``\SIb{h}{cm}``). Four +partitioning the potential transpiration rate ``T_p`` based on the fraction of the total root +length (`rootfraction` ``\SIb{}{-}``) in each soil layer. A root water uptake reduction model +is used to calculate a reduction coefficient as a function of soil water pressure, that may +reduce the maximum possible root water extraction rate. The root water uptake reduction model +is based on the concept proposed by Feddes et al. (1978). This concept defines a reduction +coefficient ``\SIb{\alpha}{-}`` as a function of soil water pressure (``\SIb{h}{cm}``). Four different levels of ``h`` are defined: `h1`, `h2`, `h3` and `h4`. `h1` represents anoxic -moisture conditions, `h2` represents field capacity, `h3` represents the point of critical -soil moisture content (onset of drought stress), and `h4` represents the wilting point. The -value of `h3` is a function of the potential transpiration rate, between ``\SIb{1}{mm d^{-1}}`` and ``\SIb{5}{mm d^{-1}}``. If ``T_p \le \SIb{1}{mm d^{-1}}``, `h3` is set equal to `h3_low` (input model -parameter). If ``T_p \ge \SIb{5}{mm d^{-1}}``, `h3` is set equal to `h3_high` (input model -parameter). For ``T_p`` values between ``\SIb{1}{mm d^{-1}}`` and ``\SIb{5}{mm d^{-1}}``, the value of `h3` is linearly -related to ``T_p`` (between `h3_low` and `h3_high`). Besides model parameters `h3_high` and -`h3_low`, the critical pressure heads `h1`, `h2` and `h4` can be defined as input to the -model. +moisture conditions, `h2` represents field capacity, `h3` represents the point of critical soil +moisture content (onset of drought stress), and `h4` represents the wilting point. The value of +`h3` is a function of the potential transpiration rate, between ``\SIb{1}{mm d^{-1}}`` and +``\SIb{5}{mm d^{-1}}``. If ``T_p \le \SIb{1}{mm d^{-1}}``, `h3` is set equal to `h3_low` (input +model parameter). If ``T_p \ge \SIb{5}{mm d^{-1}}``, `h3` is set equal to `h3_high` (input +model parameter). For ``T_p`` values between ``\SIb{1}{mm d^{-1}}`` and ``\SIb{5}{mm d^{-1}}``, +the value of `h3` is linearly related to ``T_p`` (between `h3_low` and `h3_high`). Besides +model parameters `h3_high` and `h3_low`, the critical pressure heads `h1`, `h2` and `h4` can be +defined as input to the model. -The current soil water pressure is determined following the concept defined by Brooks and -Corey (1964): +The current soil water pressure is determined following the concept defined by Brooks and Corey +(1964): $$ \frac{\theta-\theta_r}{\theta_s-\theta_r} = \min\left\{1, \left(\frac{h_b}{h}\right)^\lambda\right\} @@ -213,12 +212,11 @@ $\theta$, $\theta_s$, $\theta_r$ and $\lambda$ as previously defined. Whenever the current soil water pressure drops below `h4`, the root water uptake is set to zero. The root water uptake is at ideal conditions whenever the soil water pressure is above -`h3`, with a linear transition between `h3` and `h4`. The assumption that very wet -conditions do not affect root water uptake too much is probably generally applicable to -natural vegetation. For crops this assumption is not valid and in this case root water -uptake above `h1` should be set to zero (oxygen deficit) and between `h1` and `h2` root -water uptake is limited. This is possible by setting the input model parameter `alpha_h1` at -0 (default is 1). +`h3`, with a linear transition between `h3` and `h4`. The assumption that very wet conditions +do not affect root water uptake too much is probably generally applicable to natural +vegetation. For crops this assumption is not valid and in this case root water uptake above +`h1` should be set to zero (oxygen deficit) and between `h1` and `h2` root water uptake is +limited. This is possible by setting the input model parameter `alpha_h1` at 0 (default is 1). ![Root water uptake reduction coefficient as a function of soil water pressure](../../images/soil_rootwateruptake.png) @@ -228,11 +226,11 @@ water uptake is limited. This is possible by setting the input model parameter ` --> The maximum allowed root water extraction from each soil layer in the unsaturated zone is -determined based on the fraction of each soil layer in the unsaturated zone that is above -the rooting depth (`availcap`) and the unsaturated storage `usld`, see conceptual figure -below. This is implemented using the following code (`i` refers to the index of the vector -that contains all active cells within the spatial model domain and `k` refers to the soil -layer (from top to bottom) in the unsaturated zone): +determined based on the fraction of each soil layer in the unsaturated zone that is above the +rooting depth (`availcap`) and the unsaturated storage `usld`, see conceptual figure below. +This is implemented using the following code (`i` refers to the index of the vector that +contains all active cells within the spatial model domain and `k` refers to the soil layer +(from top to bottom) in the unsaturated zone): ```julia # availcap is fraction of soil layer containing roots @@ -255,9 +253,8 @@ layer (from top to bottom) in the unsaturated zone): ::: {.callout-note} -When `whole_ust_available` is set to true in the TOML file, almost the complete -unsaturated storage (99%) is available for transpiration, independent of the -`rootingdepth`. +When `whole_ust_available` is set to true in the TOML file, almost the complete unsaturated +storage (99%) is available for transpiration, independent of the `rootingdepth`. ```toml [model] @@ -267,19 +264,19 @@ whole_ust_available = true The computation of transpiration from the saturated store depends on the water table depth, rooting depth, the reduction coefficient $\alpha$, the fraction of wet roots and the -`rootfraction` below the water table. The fraction of wet roots is determined using a -sigmoid fuction (see figure below). The parameter `rootdistpar` defines the sharpness of the -transition between fully wet and fully dry roots. If the water table depth is equal to or -lower than the rooting depth, the remaining potential transpiration is used based on the -potential transpiration and actual transpiration in the unsaturated zone. The remaining -potential transpiration is multiplied by the wet roots fraction and the reduction -coefficient (and limited by the available water in saturated zone) to get the transpiration -from the saturated part of the soil. If the water table depth intersects the rooting depth, -the potential transpiration is multiplied by the remaining `rootfraction` (below the water -table), wet roots fraction and the reduction coefficient (and limited by the available water -in saturated zone) to get the transpiration from the saturated part of the soil. This is -implemented using the following code (`i` refers to the index of the vector that contains -all active cells within the spatial model domain): +`rootfraction` below the water table. The fraction of wet roots is determined using a sigmoid +fuction (see figure below). The parameter `rootdistpar` defines the sharpness of the transition +between fully wet and fully dry roots. If the water table depth is equal to or lower than the +rooting depth, the remaining potential transpiration is used based on the potential +transpiration and actual transpiration in the unsaturated zone. The remaining potential +transpiration is multiplied by the wet roots fraction and the reduction coefficient (and +limited by the available water in saturated zone) to get the transpiration from the saturated +part of the soil. If the water table depth intersects the rooting depth, the potential +transpiration is multiplied by the remaining `rootfraction` (below the water table), wet roots +fraction and the reduction coefficient (and limited by the available water in saturated zone) +to get the transpiration from the saturated part of the soil. This is implemented using the +following code (`i` refers to the index of the vector that contains all active cells within the +spatial model domain): ```julia # transpiration from saturated store @@ -324,11 +321,11 @@ glacier = true ### The SBM soil water accounting scheme -A detailed description of the `Topog_SBM` model has been given by Vertessy (1999). Briefly: -the soil is considered as a bucket with a certain depth ($\SIb{z_t}{mm}$), divided into a -saturated store ($\SIb{S}{mm}$) and an unsaturated store ($\SIb{U}{mm}$). The top of the $S$ -store forms a pseudo-water table at depth $\SIb{z_{i}}{mm}$ such that the value of $S$ at any -time is given by: +A detailed description of the `Topog_SBM` model has been given by Vertessy (1999). Briefly: the +soil is considered as a bucket with a certain depth ($\SIb{z_t}{mm}$), divided into a saturated +store ($\SIb{S}{mm}$) and an unsaturated store ($\SIb{U}{mm}$). The top of the $S$ store forms +a pseudo-water table at depth $\SIb{z_{i}}{mm}$ such that the value of $S$ at any time is given +by: $$ S=(z_t-z_i)(\theta_s-\theta_r) @@ -351,22 +348,23 @@ $$ S_d=(\theta_s-\theta_r)z_t-S $$ -All infiltrating water that enters the $U$ store first. The unsaturated layer can be -split-up in different layers, by providing the thickness $\SIb{}{mm}$ of the layers in the TOML -file. The following example specifies three layers (from top to bottom) of 100, 300 and 800 -mm: +All infiltrating water that enters the $U$ store first. The unsaturated layer can be split-up +in different layers, by providing the thickness $\SIb{}{mm}$ of the layers in the TOML file. +The following example specifies three layers (from top to bottom) of 100, 300 and 800 mm: ```toml [model] thicknesslayers = [100, 300, 800] ``` -The code checks for each grid cell the specified layers against the `soilthickness` $\SIb{}{mm}$, -and adds or removes (partly) layer(s) based on the `soilthickness`. +The code checks for each grid cell the specified layers against the `soilthickness` +$\SIb{}{mm}$, and adds or removes (partly) layer(s) based on the `soilthickness`. -Assuming a unit head gradient, the transfer of water ($\SIb{\mathrm{st}}{mm t^{-1}}$) from a $\SIb{U}{mm}$ store layer is controlled by the saturated hydraulic conductivity $\SIb{\subtext{K}{sat}}{mm t^{-1}}$ at depth $\SIb{z}{mm}$ (bottom layer) or $\SIb{z_i}{mm}$, the effective saturation -degree of the layer, and a Brooks-Corey power coefficient (parameter $c$) based on the -pore size distribution index $\lambda$ (Brooks and Corey, 1964): +Assuming a unit head gradient, the transfer of water ($\SIb{\mathrm{st}}{mm t^{-1}}$) from a +$\SIb{U}{mm}$ store layer is controlled by the saturated hydraulic conductivity +$\SIb{\subtext{K}{sat}}{mm t^{-1}}$ at depth $\SIb{z}{mm}$ (bottom layer) or $\SIb{z_i}{mm}$, +the effective saturation degree of the layer, and a Brooks-Corey power coefficient (parameter +$c$) based on the pore size distribution index $\lambda$ (Brooks and Corey, 1964): $$ \begin{gathered} @@ -385,14 +383,17 @@ original `Topog\_SBM` vertical transfer formulation, by specifying in the TOML f transfermethod = true ``` -The transfer of water from the $\SIb{U}{mm}$ store to the $\SIb{S}{mm}$ store ($\SIb{st}{mm t^{-1}}$) is in that case controlled by the saturated hydraulic conductivity $\SIb{\subtext{K}{sat}}{mm s^{-1}}$ at depth $\SIb{z_i}{mm}$ and the ratio between $\SIb{U}{mm}$ and $\SIb{S_d}{mm}$: +The transfer of water from the $\SIb{U}{mm}$ store to the $\SIb{S}{mm}$ store ($\SIb{st}{mm +t^{-1}}$) is in that case controlled by the saturated hydraulic conductivity +$\SIb{\subtext{K}{sat}}{mm s^{-1}}$ at depth $\SIb{z_i}{mm}$ and the ratio between +$\SIb{U}{mm}$ and $\SIb{S_d}{mm}$: $$ \mathrm{st}=\subtext{K}{sat}\frac{U_s}{S_d} $$ -Four different saturated hydraulic conductivity depth profiles (`ksat_profile`) are -available and a `ksat_profile` can be specified in the TOML file as follows: +Four different saturated hydraulic conductivity depth profiles (`ksat_profile`) are available +and a `ksat_profile` can be specified in the TOML file as follows: ```toml [input.vertical] @@ -400,21 +401,22 @@ ksat_profile = "exponential_constant" # optional, one of ("exponential", "expone ``` Soil measurements are often available for about the upper $\SI{1.5-2}{m}$ of the soil column to -estimate the saturated hydraulic conductivity, while these measurements are often lacking -for soil depths beyond $\SI{1.5-2}{m}$. These different profiles allow to extent the saturated +estimate the saturated hydraulic conductivity, while these measurements are often lacking for +soil depths beyond $\SI{1.5-2}{m}$. These different profiles allow to extent the saturated hydraulic conductivity profile based on measurements (either an exponential fit or hydraulic conductivity value per soil layer) with an exponential or constant profile. By default, with -`ksat_profile` "exponential", the saturated hydraulic conductivity $\SIb{\subtext{K}{sat}}{mm t^{-1}}$ declines with soil depth $\SIb{z}{mm}$ in the model according to: +`ksat_profile` "exponential", the saturated hydraulic conductivity $\SIb{\subtext{K}{sat}}{mm +t^{-1}}$ declines with soil depth $\SIb{z}{mm}$ in the model according to: $$ \subtext{K}{sat} = K_0 e^{-fz}, $$ -where $\SIb{K_0}{mm s^{-1}}$ is the saturated hydraulic conductivity at the soil surface -and $\SIb{f}{mm^{-1}}$ is a scaling parameter. +where $\SIb{K_0}{mm s^{-1}}$ is the saturated hydraulic conductivity at the soil surface and +$\SIb{f}{mm^{-1}}$ is a scaling parameter. -The plot below shows the relation between soil depth $z$ and saturated hydraulic -conductivity $\subtext{K}{sat}$ for different values of $f$. +The plot below shows the relation between soil depth $z$ and saturated hydraulic conductivity +$\subtext{K}{sat}$ for different values of $f$. ```{julia} # | code-fold: true @@ -422,25 +424,25 @@ using Printf using CairoMakie let - fig = Figure(resolution = (800, 400)) - ax = Axis(fig[1, 1], xlabel = L"K_\mathrm{sat}\;[\mathrm{mm/day}]", ylabel = L"-z\;[\mathrm{mm}]") + fig = Figure(resolution=(800, 400)) + ax = Axis(fig[1, 1], xlabel=L"K_\mathrm{sat}\;[\mathrm{mm/day}]", ylabel=L"-z\;[\mathrm{mm}]") - z = 0:5.0:1000 - ksat = 100.0 + z = 0:5.0:1000 + ksat = 100.0 f = 0.6 ./ collect(50:150.0:800) - for fi in f - lines!(ax, ksat .* exp.(-fi .* z), -z, label = @sprintf("%.2e", fi)) - end + for fi in f + lines!(ax, ksat .* exp.(-fi .* z), -z, label=@sprintf("%.2e", fi)) + end - Legend(fig[1, 2], ax, L"f") + Legend(fig[1, 2], ax, L"f") fig -end +end ``` -With `ksat_profile` "exponential\_constant", $\subtext{K}{sat}$ declines exponentially with soil -depth $\SIb{z}{mm}$ until $\SIb{\subtext{z}{mm}}{mm}$ below the soil surface, and stays constant at and -beyond soil depth $\subtext{z}{exp}$: +With `ksat_profile` "exponential\_constant", $\subtext{K}{sat}$ declines exponentially with +soil depth $\SIb{z}{mm}$ until $\SIb{\subtext{z}{mm}}{mm}$ below the soil surface, and stays +constant at and beyond soil depth $\subtext{z}{exp}$: $$ \subtext{K}{sat} = \begin{cases} @@ -450,38 +452,39 @@ $$ $$ It is also possible to provide a $\subtext{K}{sat}$ value per soil layer by specifying -`ksat_profile` "layered", these $\subtext{K}{sat}$ values are used directly to compute the vertical -transfer of water between soil layers and to the saturated store $S$. Finally, with the -`ksat_profile` "layered\_exponential" a $\subtext{K}{sat}$ value per soil layer is used until depth -$\subtext{z}{layered}$ below the soil surface, and beyond $\subtext{z}{layered}$ an -exponential decline of $\subtext{K}{sat}$ (of the soil layer with bottom $\subtext{z}{layered}$) -controlled by $f$ occurs. The different available `ksat_profle` options are schematized in -the figure below where the blue line represents the $\subtext{K}{sat}$ value. +`ksat_profile` "layered", these $\subtext{K}{sat}$ values are used directly to compute the +vertical transfer of water between soil layers and to the saturated store $S$. Finally, with +the `ksat_profile` "layered\_exponential" a $\subtext{K}{sat}$ value per soil layer is used +until depth $\subtext{z}{layered}$ below the soil surface, and beyond $\subtext{z}{layered}$ an +exponential decline of $\subtext{K}{sat}$ (of the soil layer with bottom +$\subtext{z}{layered}$) controlled by $f$ occurs. The different available `ksat_profle` options +are schematized in the figure below where the blue line represents the $\subtext{K}{sat}$ +value. ![Overview of available `ksat_profile` options, for a soil column with five layers](../../images/sbm_ksat_profiles.png) ### Infiltration -The water available for infiltration is taken as the rainfall including meltwater. -Infiltration is determined separately for the compacted and non-compacted areas, as these -have different infiltration capacities. Naturally, only the water that can be stored in the -soil can infiltrate. If not all water can infiltrate, this is added as excess water to the -runoff routing scheme. +The water available for infiltration is taken as the rainfall including meltwater. Infiltration +is determined separately for the compacted and non-compacted areas, as these have different +infiltration capacities. Naturally, only the water that can be stored in the soil can +infiltrate. If not all water can infiltrate, this is added as excess water to the runoff +routing scheme. The infiltrating water is split in two parts, the part that falls on compacted areas and the part that falls on non-compacted areas. The maximum amount of water that can infiltrate in these areas is calculated by taking the minimum of the maximum infiltration rate -(`infiltcapsoil` $\SIb{}{mm t^{-1}}$ for non-compacted areas and `infiltcappath` $\SIb{}{mm t^{-1}}$ -for compacted areas) and the amount of water available for infiltration `avail_forinfilt` -$\SIb{}{mm t^{-1}}$. The water that can actually infiltrate `infiltsoilpath` $\SIb{}{mm t^{-1}}$ is -calculated by taking the minimum of the total maximum infiltration rate (compacted and -non-compacted areas) and the remaining storage capacity. - -Infiltration excess occurs when the infiltration capacity is smaller then the throughfall -and stemflow rate. This amount of water (`infiltexcess` $\SIb{}{mm t^{-1}}$) becomes overland -flow (infiltration excess overland flow). Saturation excess occurs when the (upper) soil -becomes saturated and water cannot infiltrate anymore. This amount of water `excesswater` -$\SIb{}{mm t^{-1}}$ becomes overland flow (saturation excess overland flow). +(`infiltcapsoil` $\SIb{}{mm t^{-1}}$ for non-compacted areas and `infiltcappath` $\SIb{}{mm +t^{-1}}$ for compacted areas) and the amount of water available for infiltration +`avail_forinfilt` $\SIb{}{mm t^{-1}}$. The water that can actually infiltrate `infiltsoilpath` +$\SIb{}{mm t^{-1}}$ is calculated by taking the minimum of the total maximum infiltration rate +(compacted and non-compacted areas) and the remaining storage capacity. + +Infiltration excess occurs when the infiltration capacity is smaller then the throughfall and +stemflow rate. This amount of water (`infiltexcess` $\SIb{}{mm t^{-1}}$) becomes overland flow +(infiltration excess overland flow). Saturation excess occurs when the (upper) soil becomes +saturated and water cannot infiltrate anymore. This amount of water `excesswater` $\SIb{}{mm +t^{-1}}$ becomes overland flow (saturation excess overland flow). #### Infiltration in frozen soils @@ -492,18 +495,18 @@ the soil is frozen. The soil temperature is calculated based on the soil tempera previous timestep, and the temperature difference between air and soil temperature weighted with a factor (`w_soil`, which defaults to 0.1125). -The near surface soil temperature is modelled using a simple equation (Wigmosta et al., -2009): +The near surface soil temperature is modelled using a simple equation (Wigmosta et al., 2009): $$ T_s^t = T_s^{t-1} + w (T_a - T_s^{t-1}) $$ -where $\SIb{T_s^{t}}{\degree C}$ is the near-surface soil temperature at time $t$, $\SIb{T_a}{\degree C}$ is air temperature and $\SIb{w}{-}$ is a weighting coefficient determined -through calibration (default is 0.1125 for daily timesteps). +where $\SIb{T_s^{t}}{\degree C}$ is the near-surface soil temperature at time $t$, +$\SIb{T_a}{\degree C}$ is air temperature and $\SIb{w}{-}$ is a weighting coefficient +determined through calibration (default is 0.1125 for daily timesteps). -A reduction factor (`cf_soil` $\SIb{}{-}$, default is 0.038) is applied to the maximum infiltration -rate (`infiltcapsoil` and `infiltcappath`), when the following model settings are specified -in the TOML file: +A reduction factor (`cf_soil` $\SIb{}{-}$, default is 0.038) is applied to the maximum +infiltration rate (`infiltcapsoil` and `infiltcappath`), when the following model settings are +specified in the TOML file: ```toml [model] @@ -511,8 +514,8 @@ soilinfreduction = true snow = true ``` -If `soilinfreduction` is set to `false`, water is allowed to infiltrate the soil, even if -the soil is frozen. +If `soilinfreduction` is set to `false`, water is allowed to infiltrate the soil, even if the +soil is frozen. A S-curve (see plot below) is used to make a smooth transition (a c-factor ($c$) of 8.0 is used): @@ -534,8 +537,8 @@ $$ ### Capillary rise The actual capillary rise `actcapflux` $\SIb{}{mm t^{-1}}$ is determined using the following -approach: first the saturated hydraulic conductivity `ksat` $\SIb{}{mm t^{-1}}$ is determined at -the water table $z_i$; next a potential capillary rise `maxcapflux` $\SIb{}{mm t^{-1}}$ is +approach: first the saturated hydraulic conductivity `ksat` $\SIb{}{mm t^{-1}}$ is determined +at the water table $z_i$; next a potential capillary rise `maxcapflux` $\SIb{}{mm t^{-1}}$ is determined from the minimum of `ksat`, actual transpiration `actevapustore` $\SIb{}{mm t^{-1}}$ taken from the $U$ store, available water in the $S$ store (`satwaterdepth` $\SIb{}{mm}$) and the deficit of the $U$ store (`ustorecapacity` $\SIb{}{mm}$), as shown by the following code @@ -546,8 +549,9 @@ maxcapflux = max(0.0, min(ksat, actevapustore, ustorecapacity, satwaterdepth)) ``` Then the potential rise `maxcapflux` is scaled using the water table depth `zi`, a maximum -water depth `cap_hmax` $\SIb{}{mm}$ beyond which capillary rise ceases and a coefficient `cap_n` $\SIb{}{-}$, as follows in the code block below (`i` refers to the index of the vector that contains -all active cells within the spatial model domain): +water depth `cap_hmax` $\SIb{}{mm}$ beyond which capillary rise ceases and a coefficient +`cap_n` $\SIb{}{-}$, as follows in the code block below (`i` refers to the index of the vector +that contains all active cells within the spatial model domain): ```julia if sbm.zi[i] > rootingdepth @@ -561,25 +565,24 @@ else end ``` -If the roots reach the water table (`rootingdepth` $\ge$ `sbm.zi`), `capflux` is set to -zero. +If the roots reach the water table (`rootingdepth` $\ge$ `sbm.zi`), `capflux` is set to zero. Finally, the capillary rise `capflux` is limited by the unsaturated store deficit (one or -multiple layers), calculated as follows in the code block below (`i` refers to the index of -the vector that contains all active cells within the spatial model domain, and `k` refers to -the layer position): +multiple layers), calculated as follows in the code block below (`i` refers to the index of the +vector that contains all active cells within the spatial model domain, and `k` refers to the +layer position): ```julia usl[k] * (sbm.theta_s[i] - sbm.theta_r[i]) - usld[k] ``` -where `usl` $\SIb{}{mm}$ is the unsaturated layer thickness, `usld` is the `ustorelayerdepth` $\SIb{}{mm}$ -(amount of water in the unsaturated layer), and $\theta_s$ and $\theta_r$ as +where `usl` $\SIb{}{mm}$ is the unsaturated layer thickness, `usld` is the `ustorelayerdepth` +$\SIb{}{mm}$ (amount of water in the unsaturated layer), and $\theta_s$ and $\theta_r$ as previously defined. -The calculation of the actual capillary rise `actcapflux` is as follows in the code block -below (`i` refers to the index of the vector that contains all active cells within the -spatial model domain, and `k` refers to the layer position): +The calculation of the actual capillary rise `actcapflux` is as follows in the code block below +(`i` refers to the index of the vector that contains all active cells within the spatial model +domain, and `k` refers to the layer position): ```julia actcapflux = 0.0 @@ -593,14 +596,14 @@ for k = n_usl:-1:1 end ``` -In case of multiple unsaturated layers (`n_usl` $>$ 1), the calculation of the actual -capillary rise starts at the lowest unsaturated layer while keeping track of the remaining -capillary rise `netcapflux` $\SIb{}{mm t^{-1}}$. +In case of multiple unsaturated layers (`n_usl` $>$ 1), the calculation of the actual capillary +rise starts at the lowest unsaturated layer while keeping track of the remaining capillary rise +`netcapflux` $\SIb{}{mm t^{-1}}$. ### Leakage -If the `maxleakage` $\SIb{}{mm day^{-1}}$ input model parameter is set > 0, water is lost from the -saturated zone and runs out of the model. +If the `maxleakage` $\SIb{}{mm day^{-1}}$ input model parameter is set > 0, water is lost from +the saturated zone and runs out of the model. ## Open water @@ -609,9 +612,9 @@ fractions of river and lakes of each grid cell. The amount of evaporation from o assumed to be equal to potential evaporation (if sufficient water is available). ## Non-irrigation -Non-irrigation water demand and allocation computations are supported for the sectors -domestic, industry and livestock. These computations can be enabled by specifying the -following in the TOML file: +Non-irrigation water demand and allocation computations are supported for the sectors domestic, +industry and livestock. These computations can be enabled by specifying the following in the +TOML file: ```toml [model.water_demand] @@ -625,16 +628,16 @@ net demand ($d_\mathrm{net}$ $\SIb{}{mm t^{-1}}$) are provided to the model (inp cyclic or forcing data). Gross demand represents the total demand and hence the total abstraction from surface water or groundwater when sufficient water is available. Net demand represents water consumption. The portion of total abstracted water that is not consumed is -returned as surface water. The return flow fraction ($f_\mathrm{return}$ [-]) is -calculated as follows: +returned as surface water. The return flow fraction ($f_\mathrm{return}$ [-]) is calculated as +follows: $$ \subtext{f}{return} = 1.0 - \frac{\subtext{d}{net}}{\subtext{d}{gross}}, $$ -and used to calculate the return flow rate (water abstracted from surface water or -groundwater but not consumed). For grid cells containing a river the return flow is directly -returned to the river routing component, otherwise the return flow is returned to the -overland flow routing component. +and used to calculate the return flow rate (water abstracted from surface water or groundwater +but not consumed). For grid cells containing a river the return flow is directly returned to +the river routing component, otherwise the return flow is returned to the overland flow routing +component. ## Non-paddy irrigation Non-paddy (other crops than flooded rice) water demand and allocation computations are @@ -651,48 +654,50 @@ $$ (\subtext{U}{field} - \subtext{U}{a}) \ge (\subtext{U}{field} - \subtext{U}{h3}) $$ where $\SIb{\subtext{U}{field}}{mm}$ is the unsaturated store in the root zone at field -capacity (defined at a soil water pressure head of $\SI{-100}{cm}$), $\SIb{\subtext{U}{a}}{mm}$ is the -actual unsaturated store in the root zone and $\SIb{\subtext{U}{h3}}{mm}$ is the unsaturated -store in the root zone at the critical soil water pressure head `h3`, below this pressure -head reduction of root water uptake starts due to drought stress. The net irrigation demand -$\SIb{}{mm t^{-1}}$ is the irrigation rate that brings the root zone back to field capacity, -limited by the soil infiltration capacity $\SIb{}{mm t^{-1}}$, assuming that farmers do not -apply an irrigation rate higher than the soil infiltration capacity. To account for limited -irrigation efficiency the net irrigation demand is divided by the irrigation efficiency for -non-paddy crops (`irrigation_efficiency` $\SIb{}{-}$, default is $1.0$), resulting in gross irrigation -demand $\SIb{}{mm t^{-1}}$. Finally, the gross irrigation demand is limited by the maximum -irrigation rate (`maximum_irrigation_rate` $\SIb{}{mm t^{-1}}$, default is $\SI{25}{mm\;day-1}$). If -the maximum irrigation rate is applied, irrigation continues at subsequent time steps until -field capacity is reached. Irrigation is added to the `SBM` variable `avail_forinfilt` $\SIb{}{mm t^{-1}}$, the amount of water available for infiltration. +capacity (defined at a soil water pressure head of $\SI{-100}{cm}$), $\SIb{\subtext{U}{a}}{mm}$ +is the actual unsaturated store in the root zone and $\SIb{\subtext{U}{h3}}{mm}$ is the +unsaturated store in the root zone at the critical soil water pressure head `h3`, below this +pressure head reduction of root water uptake starts due to drought stress. The net irrigation +demand $\SIb{}{mm t^{-1}}$ is the irrigation rate that brings the root zone back to field +capacity, limited by the soil infiltration capacity $\SIb{}{mm t^{-1}}$, assuming that farmers +do not apply an irrigation rate higher than the soil infiltration capacity. To account for +limited irrigation efficiency the net irrigation demand is divided by the irrigation efficiency +for non-paddy crops (`irrigation_efficiency` $\SIb{}{-}$, default is $1.0$), resulting in gross +irrigation demand $\SIb{}{mm t^{-1}}$. Finally, the gross irrigation demand is limited by the +maximum irrigation rate (`maximum_irrigation_rate` $\SIb{}{mm t^{-1}}$, default is +$\SI{25}{mm\;day-1}$). If the maximum irrigation rate is applied, irrigation continues at +subsequent time steps until field capacity is reached. Irrigation is added to the `SBM` +variable `avail_forinfilt` $\SIb{}{mm t^{-1}}$, the amount of water available for infiltration. ## Paddy irrigation -Paddy (flooded rice) water demand and allocation computations are supported. These -computations can be enabled by specifying the following in the TOML file: +Paddy (flooded rice) water demand and allocation computations are supported. These computations +can be enabled by specifying the following in the TOML file: ```toml [model.water_demand] paddy = true ``` Irrigation is applied during the growing season (when input parameter `irrigation_trigger` -$\SIb{}{-}$ is `true` (or `on`)) and when the paddy water depth `h` $\SIb{}{mm}$ reaches below the minimum -water depth `h_min` $\SIb{}{mm}$ (see also the figure below). The net irrigation demand [mm -t$^{-1}$] is the irrigation rate required to reach the optimal paddy water depth `h_opt` -$\SIb{}{mm}$, an approach similar to Xie and Cui (2011). To account for limited irrigation -efficiency the net irrigation demand is divided by the irrigation efficiency for paddy -fields (`irrigation_efficiency` $\SIb{}{-}$, default is 1.0), resulting in gross irrigation demand -$\SIb{}{mm t^{-1}}$. Finally, the gross irrigation demand is limited by the maximum irrigation -rate (`maximum_irrigation_rate` $\SIb{}{mm t^{-1}}$, default is $\SIb{25}{mm d^{-1}}$). If the -maximum irrigation rate is applied, irrigation continues at subsequent time steps until the -optimal paddy water depth `h_opt` is reached. Irrigation is added to the `SBM` variable -`avail_forinfilt` $\SIb{}{mm t^{-1}}$, the amount of water available for infiltration. When the -paddy water depth `h` exceeds `h_max` $\SIb{}{mm}$ runoff occurs, and this amount is added to the -runoff routing scheme for overland flow. The figure below shows a typical vertical soil -profile of a puddled rice soil with a muddy layer of about 15 cm (in this case represented +$\SIb{}{-}$ is `true` (or `on`)) and when the paddy water depth `h` $\SIb{}{mm}$ reaches below +the minimum water depth `h_min` $\SIb{}{mm}$ (see also the figure below). The net irrigation +demand [mm t$^{-1}$] is the irrigation rate required to reach the optimal paddy water depth +`h_opt` $\SIb{}{mm}$, an approach similar to Xie and Cui (2011). To account for limited +irrigation efficiency the net irrigation demand is divided by the irrigation efficiency for +paddy fields (`irrigation_efficiency` $\SIb{}{-}$, default is 1.0), resulting in gross +irrigation demand $\SIb{}{mm t^{-1}}$. Finally, the gross irrigation demand is limited by the +maximum irrigation rate (`maximum_irrigation_rate` $\SIb{}{mm t^{-1}}$, default is $\SIb{25}{mm +d^{-1}}$). If the maximum irrigation rate is applied, irrigation continues at subsequent time +steps until the optimal paddy water depth `h_opt` is reached. Irrigation is added to the `SBM` +variable `avail_forinfilt` $\SIb{}{mm t^{-1}}$, the amount of water available for infiltration. +When the paddy water depth `h` exceeds `h_max` $\SIb{}{mm}$ runoff occurs, and this amount is +added to the runoff routing scheme for overland flow. The figure below shows a typical vertical +soil profile of a puddled rice soil with a muddy layer of about 15 cm (in this case represented by two soil layers of 5 cm and 10 cm thickness), a plow soil layer of 5 cm with relative low permeability (vertical hydraulic conductivity $k_v$ of about $\SI{5}{mm d^{-1}}$), and a -non-puddled soil below the plow soil layer. The low vertical hydraulic conductivity of the -plow soil layer can be realized by making use of the parameter `kvfrac` $\SIb{}{-}$, a -multiplication factor applied to the vertical hydraulic conductivity at soil depth $\SIb{z}{mm}$. +non-puddled soil below the plow soil layer. The low vertical hydraulic conductivity of the plow +soil layer can be realized by making use of the parameter `kvfrac` $\SIb{}{-}$, a +multiplication factor applied to the vertical hydraulic conductivity at soil depth +$\SIb{z}{mm}$. ![paddy_profile](../../images/paddy_profile.png) @@ -701,24 +706,25 @@ multiplication factor applied to the vertical hydraulic conductivity at soil dep ## Water withdrawal and allocation For the water withdrawal the total gross demand is computed (sum over the irrigation and non-irrigation water demand sectors), in case sufficient water is available the water -withdrawal is equal to the total gross demand. In case of insufficient water availability, -the water withdrawal is scaled down to the available water, and allocation is then -proportional to the gross demand per sector (industry, domestic, livestock and irrigation). -Water can be abstracted from the following sources: +withdrawal is equal to the total gross demand. In case of insufficient water availability, the +water withdrawal is scaled down to the available water, and allocation is then proportional to +the gross demand per sector (industry, domestic, livestock and irrigation). Water can be +abstracted from the following sources: - surface water from rivers (max 80% of total available water) - reservoirs and lakes (max 98% of total available water) - groundwater (max 75% of total available water) -The model parameter `frac_sw_used` (fraction surface water used, default is 1.0) determines -how much water is supplied by available surface water and groundwater. +The model parameter `frac_sw_used` (fraction surface water used, default is 1.0) determines how +much water is supplied by available surface water and groundwater. ### Local First, surface water abstraction (excluding reservoir and lake locations) is computed to -satisfy local (same grid cell) water demand. The available surface water volume is limited -by a fixed scaling factor of $0.8$ to prevent rivers from completely drying out. It is assumed -that the water demand cannot be satisfied completely from local surface water and -groundwater. The next step is to satisfy the remaining water demand for allocation `areas` $\SIb{}{-}$, described in the next sub-section. +satisfy local (same grid cell) water demand. The available surface water volume is limited by a +fixed scaling factor of $0.8$ to prevent rivers from completely drying out. It is assumed that +the water demand cannot be satisfied completely from local surface water and groundwater. The +next step is to satisfy the remaining water demand for allocation `areas` $\SIb{}{-}$, +described in the next sub-section. ### Allocation areas For allocation areas the water demand $\SIb{\subtext{V}{sw, demand}}{m^3}$ and availability @@ -728,14 +734,15 @@ limited by a fixed scaling factor of $0.98$), and the total surface water abstra $$ \subtext{V}{sw, abstraction} = \min (\subtext{V}{sw, demand}, \subtext{V}{sw, availabilty}) $$ -The fraction of available surface water that can be abstracted $\SIb{\subtext{f}{sw,abstraction}}{-}$ at the allocation area level is then: +The fraction of available surface water that can be abstracted +$\SIb{\subtext{f}{sw,abstraction}}{-}$ at the allocation area level is then: $$ \subtext{f}{sw, abstraction} = \frac{\subtext{V}{sw, abstraction}}{\subtext{V}{sw, available}} $$ -This fraction is applied to the remaining available surface water of each river cell -(including lake and reservoir locations) to compute surface water abstraction at each river -cell and to update the local surface water abstraction. +This fraction is applied to the remaining available surface water of each river cell (including +lake and reservoir locations) to compute surface water abstraction at each river cell and to +update the local surface water abstraction. The fraction of water demand that can be satisfied by available surface water $\SIb{\subtext{f}{sw, allocation}}{-}$ at the allocation area level is then: @@ -743,13 +750,14 @@ $\SIb{\subtext{f}{sw, allocation}}{-}$ at the allocation area level is then: $$ \subtext{f}{sw, allocation} = \frac{\subtext{V}{sw, abstraction}}{\subtext{V}{sw, demand}} $$ -This fraction is applied to the remaining surface water demand of each land cell to compute -the allocated surface water to each land cell. +This fraction is applied to the remaining surface water demand of each land cell to compute the +allocated surface water to each land cell. Then groundwater abstraction is computed to satisfy the remaining local water demand, where groundwater abstraction is limited by a fixed scaling factor of $0.75$ applied to the -groundwater volume. Finally, for allocation `areas` the water demand $\SIb{\subtext{V}{gw,demand}}{m^3}$ and availability $\SIb{\subtext{V}{gw, availabilty}}{m^3}$ are summed, -and the total groundwater abstraction is then: +groundwater volume. Finally, for allocation `areas` the water demand +$\SIb{\subtext{V}{gw,demand}}{m^3}$ and availability $\SIb{\subtext{V}{gw, availabilty}}{m^3}$ +are summed, and the total groundwater abstraction is then: $$ \subtext{V}{gw, abstraction} = \min(\subtext{V}{gw, demand}, \subtext{V}{gw, availabilty}) @@ -763,29 +771,30 @@ $$ This fraction is applied to the remaining available groundwater of each land cell to compute groundwater abstraction and to update the local groundwater abstraction. -The fraction of water demand that can be satisfied by available groundwater $\SIb{\subtext{f}{gw,allocation}}{-}$ at the allocation area level is then: +The fraction of water demand that can be satisfied by available groundwater +$\SIb{\subtext{f}{gw,allocation}}{-}$ at the allocation area level is then: $$ \subtext{f}{gw, allocation} = \frac{\subtext{V}{gw, abstraction}}{\subtext{V}{gw, demand}} $$ -This fraction is applied to the remaining groundwater demand of each land cell to compute -the allocated groundwater to each land cell. +This fraction is applied to the remaining groundwater demand of each land cell to compute the +allocated groundwater to each land cell. ### Abstractions -Groundwater abstraction is implemented by subtracting this amount from the `recharge` -variable of the lateral subsurface flow component (kinematic wave) or the recharge `rate` of -the groundwater flow module. Surface water `abstraction` $\SIb{}{m^3 s^{-1}}$ is divided by -the flow length `dl` $\SIb{}{m}$ and subtracted from the lateral inflow of kinematic wave routing +Groundwater abstraction is implemented by subtracting this amount from the `recharge` variable +of the lateral subsurface flow component (kinematic wave) or the recharge `rate` of the +groundwater flow module. Surface water `abstraction` $\SIb{}{m^3 s^{-1}}$ is divided by the +flow length `dl` $\SIb{}{m}$ and subtracted from the lateral inflow of kinematic wave routing scheme for river flow. For the local inertial routing scheme (river and optional floodplain routing), the surface water `abstraction` $\SIb{}{m^3 s^{-1}}$ is subtracted as part of the -continuity equation of the local inertial model. For reservoir and lake locations surface -water is abstracted (`act_surfacewater_abst_vol` $\SIb{}{m^3 s^{-1}}$) from the reservoir -`volume` $\SIb{}{m^3}$ and lake `storage` $\SIb{}{m^3}$ respectively, with a subsequent update of -the lake `waterlevel` $\SIb{}{m}$. +continuity equation of the local inertial model. For reservoir and lake locations surface water +is abstracted (`act_surfacewater_abst_vol` $\SIb{}{m^3 s^{-1}}$) from the reservoir `volume` +$\SIb{}{m^3}$ and lake `storage` $\SIb{}{m^3}$ respectively, with a subsequent update of the +lake `waterlevel` $\SIb{}{m}$. ## References -+ Brooks, R. H., and Corey, A. T., 1964, Hydraulic properties of porous media, Hydrology - Papers 3, Colorado State University, Fort Collins, 27 p. ++ Brooks, R. H., and Corey, A. T., 1964, Hydraulic properties of porous media, Hydrology Papers + 3, Colorado State University, Fort Collins, 27 p. + Feddes, R.A., Kowalik, P.J. and Zaradny, H., 1978, Simulation of field water use and crop yield, Pudoc, Wageningen, Simulation Monographs. + Gash, J. H. C., 1979, An analytical model of rainfall interception by forests, Q. J. Roy. @@ -800,10 +809,10 @@ the lake `waterlevel` $\SIb{}{m}$. vegetation of variable density using an adapted analytical model, Part 2, Model validation for a tropical upland mixed cropping system, J. Hydr., 247, 239–262. + Vertessy, R., and Elsenbeer, H., 1999, Distributed modeling of storm flow generation in an - amazonian rain forest catchment: effects of model parameterization, Water Resour. Res., - 35, 2173–2187. doi: 10.1029/1999WR9000511257. + amazonian rain forest catchment: effects of model parameterization, Water Resour. Res., 35, + 2173–2187. doi: 10.1029/1999WR9000511257. + Wigmosta, M. S., Lane, L. J., Tagestad, J. D., and Coleman A. M., 2009, Hydrologic and - erosion models to assess land use and management practices affecting soil erosion, J. - Hydrol. Eng., 14, 27-41. -+ Xie, X. and Cui, Y., 2011, Development and test of SWAT for modeling hydrological - processes in irrigation districts with paddy rice, J. Hydrol., 396, pp. 61-71. + erosion models to assess land use and management practices affecting soil erosion, J. Hydrol. + Eng., 14, 27-41. ++ Xie, X. and Cui, Y., 2011, Development and test of SWAT for modeling hydrological processes + in irrigation districts with paddy rice, J. Hydrol., 396, pp. 61-71. diff --git a/docs/model_docs/vertical/sediment.qmd b/docs/model_docs/vertical/sediment.qmd index db1ab7095..f1c9a1220 100644 --- a/docs/model_docs/vertical/sediment.qmd +++ b/docs/model_docs/vertical/sediment.qmd @@ -2,14 +2,14 @@ title: Sediment --- -Over the land, soil erosion, also called soil loss, is closely linked to the water cycle. -The main processes governing sediment generation are splash erosion from rain droplets, and -sheet and rill erosion from the shear stress caused by overland flow. The intensity of soil -erosion by rain or flow depends on the land and soil characteristics such as slope, land use -or soil type. Once soil is eroded, the detached particles can be transported downslope by -overland flow. Along the transport pathways, soil particles can also be deposited due to a -low flow velocity, a change of topography in depressions, footslopes or valley bottoms, -and/or can be filtered and stopped by a change in vegetation such as field boundaries. +Over the land, soil erosion, also called soil loss, is closely linked to the water cycle. The +main processes governing sediment generation are splash erosion from rain droplets, and sheet +and rill erosion from the shear stress caused by overland flow. The intensity of soil erosion +by rain or flow depends on the land and soil characteristics such as slope, land use or soil +type. Once soil is eroded, the detached particles can be transported downslope by overland flow. +Along the transport pathways, soil particles can also be deposited due to a low flow velocity, +a change of topography in depressions, footslopes or valley bottoms, and/or can be filtered and +stopped by a change in vegetation such as field boundaries. The inland part of the sediment gathers these different processes, separated in a vertical structure for the soil loss and lateral structure for the transport in overland flow. @@ -18,9 +18,9 @@ structure for the soil loss and lateral structure for the transport in overland ## Soil Erosion The first process to consider in sediment dynamics is the generation of sediments by land -erosion. The main processes behind soil loss are rainfall erosion and overland flow erosion. -In order to model such processes at a fine time and space scale, physics-based models such -as ANSWERS and EUROSEM were chosen here. +erosion. The main processes behind soil loss are rainfall erosion and overland flow erosion. In +order to model such processes at a fine time and space scale, physics-based models such as +ANSWERS and EUROSEM were chosen here. The choice of rainfall erosion method is set up in the model section of the TOML: ```toml @@ -30,28 +30,28 @@ rainerosmethod = "answers" # Rainfall erosion equation: ["answers", "eurosem"] ### Rainfall erosion In wflow\_sediment, rainfall erosion can both be modelled using EUROSEM or ANSWERS equation. -The main difference between the models is that EUROSEM uses a more physics-based approach -based on the kinetic energy of the rain drops impacting the soil (Morgan et al, 1998), while -ANSWERS is more empirical and uses parameters from the USLE model (Beasley et al, 1991). - -In EUROSEM, rainfall erosion is modelled according to rainfall intensity and its kinetic -energy when it reaches the soil according to equations developed by Brandt (1990). As the -intensity of the rain kinetic energy depends on the length of the fall, rainfall intercepted -by vegetation will then be reduced compared to direct throughfall. The kinetic energy of -direct throughfall is estimated by (Morgan et al, 1998): +The main difference between the models is that EUROSEM uses a more physics-based approach based +on the kinetic energy of the rain drops impacting the soil (Morgan et al, 1998), while ANSWERS +is more empirical and uses parameters from the USLE model (Beasley et al, 1991). + +In EUROSEM, rainfall erosion is modelled according to rainfall intensity and its kinetic energy +when it reaches the soil according to equations developed by Brandt (1990). As the intensity of +the rain kinetic energy depends on the length of the fall, rainfall intercepted by vegetation +will then be reduced compared to direct throughfall. The kinetic energy of direct throughfall +is estimated by (Morgan et al, 1998): $$ \subtext{\mathrm{KE}}{direct} = 8.95 + 8.44\,\log_{10}(R_i) $$ -where $\SIb{\subtext{\mathrm{KE}}{direct}}{J m^{-2} mm^{-1}}$ is the kinetic energy of direct throughfall and -$\SIb{R_i}{mm h^{-1}}$ is rainfall intensity. If the rainfall is intercepted by -vegetation and falls as leaf drainage, its kinetic energy is then reduced according to -(Brandt, 1990): +where $\SIb{\subtext{\mathrm{KE}}{direct}}{J m^{-2} mm^{-1}}$ is the kinetic energy of direct +throughfall and $\SIb{R_i}{mm h^{-1}}$ is rainfall intensity. If the rainfall is intercepted by +vegetation and falls as leaf drainage, its kinetic energy is then reduced according to (Brandt, +1990): $$ \subtext{\mathrm{KE}}{leaf} = 15.8\,\sqrt{H_p} - 5.87 $$ -where $\SIb{\subtext{\mathrm{KE}}{leaf}}{J m^{-2} mm^{-1}}$ is kinetic energy of leaf drainage and -$\SIb{H_p}{m}$ is the effective canopy height (half of plant height). Canopy height can be +where $\SIb{\subtext{\mathrm{KE}}{leaf}}{J m^{-2} mm^{-1}}$ is kinetic energy of leaf drainage +and $\SIb{H_p}{m}$ is the effective canopy height (half of plant height). Canopy height can be derived from the global map from Simard & al. (2011) or by user input depending on the land use. @@ -62,14 +62,15 @@ then: $$ D_R = k\,\mathrm{KE}\,e^{-\varphi h} $$ -where $\SIb{k}{g J^{-1}}$ is an index of the detachability of the soil, $\SIb{\mathrm{KE}}{J m^{-2}}$ is the total -rainfall kinetic energy, $\SIb{h}{m}$ is the surface runoff depth on the soil and $\varphi$ is an exponent varying between $0.9$ and $3.1$ used to reduce rainfall impact if -the soil is already covered by water. As a simplification, Torri (1987) has shown that a -value of $2.0$ for $\varphi$ is representative enough for a wide range of soil conditions. -The detachability of the soil $k$ depends on the soil texture (proportion of clay, silt -and sand content) and corresponding values are defined in EUROSEM user guide (Morgan et al, -1998). As a simplification, in `wflow_sediment`, the mean value of the detachability shown in -the table below are used. Soil texture can for example be derived from the topsoil clay and +where $\SIb{k}{g J^{-1}}$ is an index of the detachability of the soil, $\SIb{\mathrm{KE}}{J +m^{-2}}$ is the total rainfall kinetic energy, $\SIb{h}{m}$ is the surface runoff depth on the +soil and $\varphi$ is an exponent varying between $0.9$ and $3.1$ used to reduce rainfall +impact if the soil is already covered by water. As a simplification, Torri (1987) has shown +that a value of $2.0$ for $\varphi$ is representative enough for a wide range of soil +conditions. The detachability of the soil $k$ depends on the soil texture (proportion of clay, +silt and sand content) and corresponding values are defined in EUROSEM user guide (Morgan et +al, 1998). As a simplification, in `wflow_sediment`, the mean value of the detachability shown +in the table below are used. Soil texture can for example be derived from the topsoil clay and silt content from SoilGrids (Hengl et al, 2017). Table: Mean detachability of soil depending on its texture (Morgan et al, 1998). @@ -86,29 +87,27 @@ Table: Mean detachability of soil depending on its texture (Morgan et al, 1998). | Fine Sand | 3.5 | | Sand | 1.9 | -Rainfall erosion is handled differently in ANSWERS. There, the impacts of vegetation and -soil properties are handled through the USLE coefficients in the equation (Beasley et al, -1991): +Rainfall erosion is handled differently in ANSWERS. There, the impacts of vegetation and soil +properties are handled through the USLE coefficients in the equation (Beasley et al, 1991): $$ D_R = 0.108 \, \subtext{C}{USLE} \, \subtext{K}{USLE} \, A_i \, R_i^2 $$ -where $\SIb{D_R}{kg min^{-1}}$ is the soil detachment by rainfall, $\subtext{C}{USLE}$ -is the soil cover-management factor from the USLE equation, $\subtext{K}{USLE}$ is the soil +where $\SIb{D_R}{kg min^{-1}}$ is the soil detachment by rainfall, $\subtext{C}{USLE}$ is the +soil cover-management factor from the USLE equation, $\subtext{K}{USLE}$ is the soil erodibility factor from the USLE equation, $\SIb{A_i}{m^2}$ is the area of the cell and -$\SIb{R_i}{mm\;min^{-1}}$ is the rainfall intensity. There are several methods -available to estimate the $C$ and $K$ factors from the USLE. They can come from user -input maps, for example maps resulting from Panagos & al.'s recent studies for Europe -(Panagos et al, 2015) (Ballabio et al, 2016). To get an estimate of the $C$ factor -globally, the other method is to estimate $C$ values for the different land use type in -from global land cover maps (e.g. GlobCover). An example is given for the global land cover -map GlobCover, summed up in the table below, the values come from a literature study -including Panagos et al.'s review (2015), Gericke & al. (2015), Mansoor & al. (2013), Chadli -et al. (2016), de Vente et al. (2009), Borrelli et al. (2014), Yang et al. (2003) and Bosco -et al. (2015). +$\SIb{R_i}{mm\;min^{-1}}$ is the rainfall intensity. There are several methods available to +estimate the $C$ and $K$ factors from the USLE. They can come from user input maps, for example +maps resulting from Panagos & al.'s recent studies for Europe (Panagos et al, 2015) (Ballabio +et al, 2016). To get an estimate of the $C$ factor globally, the other method is to estimate +$C$ values for the different land use type in from global land cover maps (e.g. GlobCover). An +example is given for the global land cover map GlobCover, summed up in the table below, the +values come from a literature study including Panagos et al.'s review (2015), Gericke & al. +(2015), Mansoor & al. (2013), Chadli et al. (2016), de Vente et al. (2009), Borrelli et al. +(2014), Yang et al. (2003) and Bosco et al. (2015). The other methods to estimate the USLE $K$ factor are to use either topsoil composition or -topsoil geometric mean diameter. $K$ estimation from topsoil composition is estimated with -the equation developed in the EPIC model (Williams et al, 1983): +topsoil geometric mean diameter. $K$ estimation from topsoil composition is estimated with the +equation developed in the EPIC model (Williams et al, 1983): $$ \begin{gathered} \subtext{K}{USLE} = \left[ 0.2 + 0.3\exp\left(-0.0256\;\mathrm{SAN}\frac{(1-\mathrm{SIL})}{100}\right) \right] @@ -116,16 +115,16 @@ $$ \left(1-\frac{0.25\;\mathrm{OC}}{\mathrm{OC}+e^{3.72-2.95\;\mathrm{OC}}}\right)\left(1-\frac{0.75\;\mathrm{SN}}{\mathrm{SN}+e^{-5.51+22.9\;\mathrm{SN}}}\right) \end{gathered} $$ -where $\SIb{\mathrm{CLA}}{\%}$, $\SIb{\mathrm{SIL}}{\%}$, $\SIb{\mathrm{SAN}}{\%}$ are respectively the clay, silt and sand fractions of the -topsoil, $\SIb{OC}{\%}$ is the topsoil organic carbon content and $\mathrm{SN} = 1-\mathrm{SAN}/100$. -These soil parameters can be derived for example from the SoilGrids dataset. The $K$ -factor can also be estimated from the soil mean geometric diameter using the formulation -from the RUSLE guide by Renard & al. (1997): +where $\SIb{\mathrm{CLA}}{\%}$, $\SIb{\mathrm{SIL}}{\%}$, $\SIb{\mathrm{SAN}}{\%}$ are +respectively the clay, silt and sand fractions of the topsoil, $\SIb{OC}{\%}$ is the topsoil +organic carbon content and $\mathrm{SN} = 1-\mathrm{SAN}/100$. These soil parameters can be +derived for example from the SoilGrids dataset. The $K$ factor can also be estimated from the +soil mean geometric diameter using the formulation from the RUSLE guide by Renard & al. (1997): $$ \subtext{K}{USLE} = 0.0034 + 0.0405\exp\left(-\dfrac{1}{2}\left(\dfrac{\log_{10}(D_g)+1.659}{0.7101}\right)^2\right) $$ -where $D_g$ is the soil geometric mean diameter (mm) estimated from topsoil clay, silt, -sand fraction. +where $D_g$ is the soil geometric mean diameter (mm) estimated from topsoil clay, silt, sand +fraction. Table: Estimation of USLE C factor per Globcover land use type @@ -157,68 +156,68 @@ Table: Estimation of USLE C factor per Globcover land use type ### Overland flow erosion -Overland flow (or surface runoff) erosion is induced by the strength of the shear stress of -the surface water on the soil. As in rainfall erosion, the effect of the flow shear stress -can be reduced by the soil vegetation or by the soil properties. In wflow_sediment, soil -detachment by overland flow is modelled as in ANSWERS with (Beasley et al, 1991): +Overland flow (or surface runoff) erosion is induced by the strength of the shear stress of the +surface water on the soil. As in rainfall erosion, the effect of the flow shear stress can be +reduced by the soil vegetation or by the soil properties. In wflow_sediment, soil detachment by +overland flow is modelled as in ANSWERS with (Beasley et al, 1991): $$ D_G = 0.90 \, \subtext{C}{USLE} \, \subtext{K}{USLE} \, A_i \, S \, q $$ -where $\SIb{D_F}{kg\;min^{-1}}$ is soil detachment by flow, $\subtext{C}{USLE}$ and $\subtext{K}{USLE}$ -are the USLE cover and soil erodibility factors, $\SIb{A_i}{m^2}$ is the cell area, -$S$ is the slope gradient and $\SIb{q}{m^2 min^{-1}}$ is the overland flow rate per unit width. The USLE $C$ and $K$ factors can be estimated with the same methods as -for rainfall erosion and here the slope gradient is obtained from the sinus rather than the -tangent of the slope angle. +where $\SIb{D_F}{kg\;min^{-1}}$ is soil detachment by flow, $\subtext{C}{USLE}$ and +$\subtext{K}{USLE}$ are the USLE cover and soil erodibility factors, $\SIb{A_i}{m^2}$ is the +cell area, $S$ is the slope gradient and $\SIb{q}{m^2 min^{-1}}$ is the overland flow rate per +unit width. The USLE $C$ and $K$ factors can be estimated with the same methods as for rainfall +erosion and here the slope gradient is obtained from the sinus rather than the tangent of the +slope angle. ## Delivery to the river system -Once soil is detached, it can be transported by overland flow and reach the river system. -This process is described in [Sediment Flux in overland flow](@ref). +Once soil is detached, it can be transported by overland flow and reach the river system. This +process is described in [Sediment Flux in overland flow](@ref). ## References + D.B Beasley and L.F Huggins. ANSWERS - Users Manual. Technical report, EPA, 1991. -+ P. Borrelli, M. Märker, P. Panagos, and B. Schütt. Modeling soil erosion and river - sediment yield for an intermountain drainage basin of the Central Apennines, Italy. - Catena, 114:45-58, 2014. 10.1016/j.catena.2013.10.007 ++ P. Borrelli, M. Märker, P. Panagos, and B. Schütt. Modeling soil erosion and river sediment + yield for an intermountain drainage basin of the Central Apennines, Italy. Catena, 114:45-58, + 2014. 10.1016/j.catena.2013.10.007 + C. Bosco, D. De Rigo, O. Dewitte, J. Poesen, and P. Panagos. Modelling soil erosion at - European scale: Towards harmonization and reproducibility. Natural Hazards and Earth - System Sciences, 15(2):225-245, 2015. 10.5194/nhess-15-225-2015 + European scale: Towards harmonization and reproducibility. Natural Hazards and Earth System + Sciences, 15(2):225-245, 2015. 10.5194/nhess-15-225-2015 + C.J Brandt. Simulation of the size distribution and erosivity of raindrops and throughfall drops. Earth Surface Processes and Landforms, 15(8):687-698, dec 1990. -+ K. Chadli. Estimation of soil loss using RUSLE model for Sebou watershed (Morocco). - Modeling Earth Systems and Environment, 2(2):51, 2016. 10.1007/s40808-016-0105-y ++ K. Chadli. Estimation of soil loss using RUSLE model for Sebou watershed (Morocco). Modeling + Earth Systems and Environment, 2(2):51, 2016. 10.1007/s40808-016-0105-y + G R Foster. Modeling the erosion process. Hydrologic modeling of small watersheds, pages 295-380, 1982. -+ A. Gericke. Soil loss estimation and empirical relationships for sediment delivery ratios - of European river catchments. International Journal of River Basin Management, 2015. ++ A. Gericke. Soil loss estimation and empirical relationships for sediment delivery ratios of + European river catchments. International Journal of River Basin Management, 2015. 10.1080/15715124.2014.1003302 + L.D.K. Mansoor, M.D. Matlock, E.C. Cummings, and L.L. Nalley. Quantifying and mapping - multiple ecosystem services change in West Africa. Agriculture, Ecosystems and - Environment, 165:6-18, 2013. 10.1016/j.agee.2012.12.001 -+ Q Morgan, J.N Smith, R.E Govers, G Poesen, J.W.A Auerswald, K Chisci, G Torri, D Styczen, - and M E Folly. The European soil erosion model (EUROSEM): documentation and user guide. - Technical report, 1998. -+ S.L Neitsch, J.G Arnold, J.R Kiniry, and J.R Williams. SWAT Theoretical Documentation - Version 2009. Texas Water Resources Institute, pages 1-647, 2011. - 10.1016/j.scitotenv.2015.11.063 + multiple ecosystem services change in West Africa. Agriculture, Ecosystems and Environment, + 165:6-18, 2013. 10.1016/j.agee.2012.12.001 ++ Q Morgan, J.N Smith, R.E Govers, G Poesen, J.W.A Auerswald, K Chisci, G Torri, D Styczen, and + M E Folly. The European soil erosion model (EUROSEM): documentation and user guide. Technical + report, 1998. ++ S.L Neitsch, J.G Arnold, J.R Kiniry, and J.R Williams. SWAT Theoretical Documentation Version + 2009. Texas Water Resources Institute, pages 1-647, 2011. 10.1016/j.scitotenv.2015.11.063 + P. Panagos, P. Borrelli, K. Meusburger, C. Alewell, E. Lugato, and L. Montanarella. - Estimating the soil erosion cover-management factor at the European scale. Land Use - Policy, 48:38-50, 2015. 10.1016/j.landusepol.2015.05.021 -+ K Renard, Gr Foster, Ga Weesies, Dk McCool, and Dc Yoder. Predicting soil erosion by - water: a guide to conservation planning with the Revised Universal Soil Loss Equation - (RUSLE). Washington, 1997. -+ D. Torri, M. Sfalanga, and M. Del Sette. Splash detachment: Runoff depth and soil - cohesion. Catena, 14(1-3):149-155, 1987. 10.1016/S0341-8162(87)80013-9 -+ J. de Vente, J. Poesen, G. Govers, and C. Boix-Fayos. The implications of data selection - for regional erosion and sediment yield modelling. Earth Surface Processes and Landforms, + Estimating the soil erosion cover-management factor at the European scale. Land Use Policy, + 48:38-50, 2015. 10.1016/j.landusepol.2015.05.021 ++ K Renard, Gr Foster, Ga Weesies, Dk McCool, and Dc Yoder. Predicting soil erosion by water: a + guide to conservation planning with the Revised Universal Soil Loss Equation (RUSLE). + Washington, 1997. ++ D. Torri, M. Sfalanga, and M. Del Sette. Splash detachment: Runoff depth and soil cohesion. + Catena, 14(1-3):149-155, 1987. 10.1016/S0341-8162(87)80013-9 ++ J. de Vente, J. Poesen, G. Govers, and C. Boix-Fayos. The implications of data selection for + regional erosion and sediment yield modelling. Earth Surface Processes and Landforms, 34(15):1994-2007, 2009. 10.1002/esp.1884 + G. Verstraeten and J. Poesen. Estimating trap efficiency of small reservoirs and ponds: methods and implications for the assessment of sediment yield. Progress in Physical Geography, 24(2):219-251, 2000. 10.1177/030913330002400204 + O. Vigiak, A. Malago, F. Bouraoui, M. Vanmaercke, and J. Poesen. Adapting SWAT hillslope - erosion model to predict sediment concentrations and yields in large Basins. Science of - the Total Environment, 538:855-875, 2015. 10.1016/j.scitotenv.2015.08.095 -+ J.R. Williams, K.G. Renard, and P.T. Dyke. EPIC A new method for assessing erosion's - effect on soil productivity. Journal of Soil and Water Conservation, 38(5):381-383, 1983. + erosion model to predict sediment concentrations and yields in large Basins. Science of the + Total Environment, 538:855-875, 2015. 10.1016/j.scitotenv.2015.08.095 ++ J.R. Williams, K.G. Renard, and P.T. Dyke. EPIC A new method for assessing erosion's effect + on soil productivity. Journal of Soil and Water Conservation, 38(5):381-383, 1983. + D. Yang, S. Kanae, T. Oki, T. Koike, and K. Musiake. Global potential soil erosion with reference to land use and climate changes. Hydrological Processes, 17(14):2913-2928, 2003. 10.1002/hyp.1441 diff --git a/docs/model_docs/vertical/shared_processes.qmd b/docs/model_docs/vertical/shared_processes.qmd index 33740ecc2..cb7aa77a2 100644 --- a/docs/model_docs/vertical/shared_processes.qmd +++ b/docs/model_docs/vertical/shared_processes.qmd @@ -7,13 +7,13 @@ title: Shared processes ### Snow modelling If the air temperature, $T_a$, is below a user-defined threshold `tt` $\SIb{}{\degree C}$ -precipitation occurs as snowfall, whereas it occurs as rainfall if $T_a ≥ \mathrm{tt}$. A another -parameter `tti` defines how precipitation can occur partly as rain or snowfall (see the +precipitation occurs as snowfall, whereas it occurs as rainfall if $T_a ≥ \mathrm{tt}$. A +another parameter `tti` defines how precipitation can occur partly as rain or snowfall (see the figure below). If precipitation occurs as snowfall, it is added to the dry snow component within the snow pack. Otherwise it ends up in the free water reservoir, which represents the liquid water content of the snow pack. Between the two components of the snow pack, -interactions take place, either through snow melt (if temperatures are above a threshold -`tt`) or through snow refreezing (if temperatures are below threshold `tt`). +interactions take place, either through snow melt (if temperatures are above a threshold `tt`) +or through snow refreezing (if temperatures are below threshold `tt`). The respective rates of snow melt and refreezing are: @@ -25,12 +25,13 @@ $$ $$ where $Q_m$ is the rate of snow melt, $Q_r$ is the rate of snow refreezing, and -$\SIb{\subtext{\mathrm{cf}}{max}}{mm\;(\degree C)^{-1} day^{-1}}$ and $\mathrm{cf}_r$ are user defined model parameters (the melting factor and the refreezing factor respectively). +$\SIb{\subtext{\mathrm{cf}}{max}}{mm\;(\degree C)^{-1} day^{-1}}$ and $\mathrm{cf}_r$ are user +defined model parameters (the melting factor and the refreezing factor respectively). The fraction of liquid water in the snow pack is at most equal to a user defined fraction, `whc`, of the water equivalent of the dry snow content. If the liquid water concentration -exceeds `whc`, either through snow melt or incoming rainfall, the surplus water -(`rainfall`) becomes available for infiltration into the soil: +exceeds `whc`, either through snow melt or incoming rainfall, the surplus water (`rainfall`) +becomes available for infiltration into the soil: ```julia snowwater = snowwater - refreezing # free water content in snow @@ -52,21 +53,21 @@ into firn/ice (using the HBV-light model) and glacier melt (using a temperature model). The definition of glacier boundaries and initial volume is defined by two parameters. The -parameter `glacierfrac` gives the fraction of each grid cell covered by a glacier as a -number between zero and one. The state parameter `glacierstore` gives the amount of water -(in mm w.e.) within the glaciers at each grid cell. Because the glacier store -(`glacierstore`) cannot be initialized by running the model for a couple of years, a default -initial state should be supplied by adding this parameter to the input static file. The -required glacier data can be prepared from available glacier datasets. +parameter `glacierfrac` gives the fraction of each grid cell covered by a glacier as a number +between zero and one. The state parameter `glacierstore` gives the amount of water (in mm w.e.) +within the glaciers at each grid cell. Because the glacier store (`glacierstore`) cannot be +initialized by running thFe model for a couple of years, a default initial state should be +supplied by adding this parameter to the input static file. The required glacier data can be +prepared from available glacier datasets. First, a fixed fraction of the snowpack on top of the glacier is converted into ice for each -timestep and added to the `glacierstore` using the HBV-light model (Seibert et al., 2018). -This fraction `g_sifrac` typically ranges from $0.001$ to $0.006$. +timestep and added to the `glacierstore` using the HBV-light model (Seibert et al., 2018). This +fraction `g_sifrac` typically ranges from $0.001$ to $0.006$. -Then, when the snowpack on top of the glacier is almost all melted (snow cover $< \SI{10}{mm}$), -glacier melt is enabled and estimated with a degree-day model. If the air temperature, -$T_a$, is below a certain threshold `g_tt` ($\SIb{}{\degree C}$) precipitation occurs as -snowfall, whereas it occurs as rainfall if $T_a ≥$ `g_tt`. +Then, when the snowpack on top of the glacier is almost all melted (snow cover $< +\SI{10}{mm}$), glacier melt is enabled and estimated with a degree-day model. If the air +temperature, $T_a$, is below a certain threshold `g_tt` ($\SIb{}{\degree C}$) precipitation +occurs as snowfall, whereas it occurs as rainfall if $T_a ≥$ `g_tt`. With this the rate of glacier melt in mm is estimated as: @@ -74,11 +75,12 @@ $$ Q_m = \subtext{g}{cfmax}(T_a − \subtext{g}{tt})\, ; \, T_a > \subtext{g}{tt} $$ -where $Q_m$ is the rate of glacier melt and $\SIb{\subtext{g}{cfmax}}{mm (\degree C)^{-1}day^{-1}}$ is the melting factor. Parameter `g_tt` can be taken as equal to the snow `tt` parameter. -Values of the melting factor `g_cfmax` normally varies from one glacier to another and some -values are reported in the literature. `g_cfmax` can also be estimated by multiplying snow -`cfmax` by a factor between 1 and 2, to take into account the higher albedo of ice compared -to snow. +where $Q_m$ is the rate of glacier melt and $\SIb{\subtext{g}{cfmax}}{mm (\degree +C)^{-1}day^{-1}}$ is the melting factor. Parameter `g_tt` can be taken as equal to the snow +`tt` parameter. Values of the melting factor `g_cfmax` normally varies from one glacier to +another and some values are reported in the literature. `g_cfmax` can also be estimated by +multiplying snow `cfmax` by a factor between 1 and 2, to take into account the higher albedo of +ice compared to snow. ## Rainfall interception Both the Gash and Rutter models are available to estimate rainfall interception by the @@ -87,8 +89,8 @@ vegetation. The selection of an interception model depends on the simulation tim ### The analytical (Gash) model The analytical model of rainfall interception is based on Rutter's numerical model. The simplifications that introduced allow the model to be applied on a daily basis, although a -storm-based approach will yield better results in situations with more than one storm per -day. The amount of water needed to completely saturate the canopy is defined as: +storm-based approach will yield better results in situations with more than one storm per day. +The amount of water needed to completely saturate the canopy is defined as: $$ P'=\frac{-\overline{R}S}{\overline{E}_w}\log\left[1-\frac{\overline{E}_w}{\overline{R}}(1-p-p_t)^{-1}\right] @@ -96,14 +98,13 @@ $$ where $\overline{R}$ is the average precipitation intensity on a saturated canopy and $\overline{E}_{w}$ the average evaporation from the wet canopy and with the vegetation -parameters $S$, $p$ and $p_t$ as defined previously. The model uses a series of -expressions to calculate the interception loss during different phases of a storm. An -analytical integration of the total evaporation and rainfall under saturated canopy -conditions is then done for each storm to determine average values of $\overline{E}_{w}$ -and $\overline{R}$. The total evaporation from the canopy (the total interception loss) is -calculated as the sum of the components listed in the table below. Interception losses from -the stems are calculated for days with $P\geq S_{t}/p_{t}$. $p_t$ and $S_t$ are small -and neglected. +parameters $S$, $p$ and $p_t$ as defined previously. The model uses a series of expressions to +calculate the interception loss during different phases of a storm. An analytical integration +of the total evaporation and rainfall under saturated canopy conditions is then done for each +storm to determine average values of $\overline{E}_{w}$ and $\overline{R}$. The total +evaporation from the canopy (the total interception loss) is calculated as the sum of the +components listed in the table below. Interception losses from the stems are calculated for +days with $P\geq S_{t}/p_{t}$. $p_t$ and $S_t$ are small and neglected. Table: Formulation of the components of interception loss according to Gash: @@ -117,16 +118,16 @@ Table: Formulation of the components of interception loss according to Gash: | Evaporation from trunks in $m+n-q$ storms that do not fill the trunk storage | $p_t\sum_{j=1}^{m+n-q}P_{g,j}$ | In applying the analytical model, saturated conditions are assumed to occur when the hourly -rainfall exceeds a certain threshold. Often a threshold of 0.5 mm/hr is used. -$\overline{R}$ is calculated for all hours when the rainfall exceeds the threshold to give -an estimate of the mean rainfall rate onto a saturated canopy. +rainfall exceeds a certain threshold. Often a threshold of 0.5 mm/hr is used. $\overline{R}$ is +calculated for all hours when the rainfall exceeds the threshold to give an estimate of the +mean rainfall rate onto a saturated canopy. -Gash (1979) has shown that in a regression of interception loss on rainfall (on a storm -basis) the regression coefficient should equal to $\overline{E}_w/\overline{R}$. Assuming -that neither $\overline{E}_w$ nor $\overline{R}$ vary considerably in time, -$\overline{E}_w$ can be estimated in this way from $\overline{R}$ in the absence of -above-canopy climatic observations. Values derived in this way generally tend to be (much) -higher than those calculated with the penman-monteith equation. +Gash (1979) has shown that in a regression of interception loss on rainfall (on a storm basis) +the regression coefficient should equal to $\overline{E}_w/\overline{R}$. Assuming that neither +$\overline{E}_w$ nor $\overline{R}$ vary considerably in time, $\overline{E}_w$ can be +estimated in this way from $\overline{R}$ in the absence of above-canopy climatic observations. +Values derived in this way generally tend to be (much) higher than those calculated with the +penman-monteith equation. ### The modified rutter model For sub daily timesteps the interception is calculated using a simplification of the Rutter @@ -139,8 +140,8 @@ Wflow.rainfall_interception_modrut ### Interception parameters from LAI The SBM concept can determine the interception parameters from leaf area index (LAI) -climatology. In order to switch this on you must define this cyclic parameter in the TOML -file, the parameter is read from `path_static`, as follows: +climatology. In order to switch this on you must define this cyclic parameter in the TOML file, +the parameter is read from `path_static`, as follows: ```toml [input] @@ -155,9 +156,9 @@ Furthermore these additional parameters are required: + Storage woody part of vegetation (`swood` $\SIb{}{mm}$) + Extinction coefficient (`kext` $\SIb{}{-}$) -Here it is assumed that `cmax` $\SIb{}{mm}$ (leaves) (canopy storage capacity for the leaves only) -relates linearly with LAI (c.f. Van Dijk and Bruijnzeel 2001). This done via the `sl`. `sl` -can be determined through a lookup table with land cover based on literature (Pitman 1989, +Here it is assumed that `cmax` $\SIb{}{mm}$ (leaves) (canopy storage capacity for the leaves +only) relates linearly with LAI (c.f. Van Dijk and Bruijnzeel 2001). This done via the `sl`. +`sl` can be determined through a lookup table with land cover based on literature (Pitman 1989, Lui 1998). Next the `cmax` (leaves) is determined using: $$ @@ -179,5 +180,5 @@ The extinction coefficient `kext` can be related to land cover. ## References + Seibert, J., Vis, M. J. P., Kohn, I., Weiler, M., and Stahl, K., 2018, Technical note: - Representing glacier geometry changes in a semi-distributed hydrological model, Hydrol. - Earth Syst. Sci., 22, 2211–2224, https://doi.org/10.5194/hess-22-2211-2018. + Representing glacier geometry changes in a semi-distributed hydrological model, Hydrol. Earth + Syst. Sci., 22, 2211–2224, https://doi.org/10.5194/hess-22-2211-2018. diff --git a/docs/user_guide/bmi.qmd b/docs/user_guide/bmi.qmd index 732fd07fd..e38a3ff4b 100644 --- a/docs/user_guide/bmi.qmd +++ b/docs/user_guide/bmi.qmd @@ -2,7 +2,6 @@ title: Basic Model Interface (BMI) --- - ## Introduction The [Community Surface Dynamics Modeling System](https://csdms.colorado.edu/wiki/Main_Page) (CSMDS) has developed the Basic Model Interface (BMI). BMI consists of a set of standard @@ -11,10 +10,9 @@ model both easier to learn and easier to couple with other software elements. For more information see also: -CSDMS provides specifications for the languages C, C++, Fortran and Python. Wflow, written -in the [Julia programming language](https://julialang.org/), makes use of the following -[Julia specification](https://github.com/Deltares/BasicModelInterface.jl), based on BMI 2.0 -version. +CSDMS provides specifications for the languages C, C++, Fortran and Python. Wflow, written in +the [Julia programming language](https://julialang.org/), makes use of the following [Julia +specification](https://github.com/Deltares/BasicModelInterface.jl), based on BMI 2.0 version. For the BMI implementation of wflow all grids are defined as [unstructured grids](https://bmi-spec.readthedocs.io/en/latest/model_grids.html#unstructured-grids), @@ -23,10 +21,10 @@ parameters) is structured (uniform rectilinear), internally wflow works with one arrays based on the active grid cells of the 2D model domain. ## Configuration -The variables that wflow can exchange through BMI are based on the different model -components and these components should be listed under the `API` section of the TOML -configuration file of the model type. Below an example of this `API` section, that lists the -`vertical` component and different `lateral` components: +The variables that wflow can exchange through BMI are based on the different model components +and these components should be listed under the `API` section of the TOML configuration file of +the model type. Below an example of this `API` section, that lists the `vertical` component and +different `lateral` components: ```toml [API] @@ -46,14 +44,14 @@ Wflow.BMI.get_input_var_names ``` Variables with a third dimension, for example `layer` as part of the vertical `SBM` concept, -are exposed as two-dimensional grids through the wflow BMI implementation. For these -variables the index of this third dimension is required, by adding `[k]` to the variable -name (`k` refers to the index of the third dimension). For example, the variable -`vertical.vwc[1]` refers to the first soil layer of the vertical `SBM` concept. +are exposed as two-dimensional grids through the wflow BMI implementation. For these variables +the index of this third dimension is required, by adding `[k]` to the variable name (`k` refers +to the index of the third dimension). For example, the variable `vertical.vwc[1]` refers to the +first soil layer of the vertical `SBM` concept. ## Couple to a groundwater model -For the coupling of wflow\_sbm (SBM + kinematic wave) with a groundwater model (e.g. -MODFLOW) it is possible to run: +For the coupling of wflow\_sbm (SBM + kinematic wave) with a groundwater model (e.g. MODFLOW) +it is possible to run: - wflow\_sbm in parts from the BMI, and - to switch off the lateral subsurface flow component of wflow\_sbm. @@ -64,11 +62,10 @@ possible to run first the recharge part of SBM: ```julia model = BMI.update(model, run="sbm_until_recharge") ``` -and to exchange recharge and for example river waterlevels to the groundwater model. After -the groundwater model update, and the exchange of groundwater head and for example drain and -river flux to wflow\_sbm, the SBM part that mainly determines exfiltration of water from the -unsaturated store, and the kinematic wave for river - and overland flow can be run as -follows: +and to exchange recharge and for example river waterlevels to the groundwater model. After the +groundwater model update, and the exchange of groundwater head and for example drain and river +flux to wflow\_sbm, the SBM part that mainly determines exfiltration of water from the +unsaturated store, and the kinematic wave for river - and overland flow can be run as follows: ```julia model = BMI.update(model, run="sbm_after_subsurfaceflow") diff --git a/docs/user_guide/faq.qmd b/docs/user_guide/faq.qmd index ddd68db3b..cd19f0a3f 100644 --- a/docs/user_guide/faq.qmd +++ b/docs/user_guide/faq.qmd @@ -4,7 +4,9 @@ title: Frequently asked questions ### How do I easily modify input parameters? -See [this section](./toml_file.qmd#modify-parameters) on how to adjust maps, and [this section](./toml_file.qmd#fixed-forcing-values) on how to directly pass uniform values. Note that both options work for any parameter. +See [this section](./toml_file.qmd#modify-parameters) on how to adjust maps, and [this +section](./toml_file.qmd#fixed-forcing-values) on how to directly pass uniform values. Note +that both options work for any parameter. ### How do I start wflow with initial conditions from a previous run?{#sec-modify-pars} @@ -12,7 +14,8 @@ See [here](./toml_file.qmd#state-options) ### How do I add external inflows and/or abstractions? -`lateral.river.inflow`: positive for inflows, negative for abstraction. If parameter is time varying, add it to the correct section, see [below](#how-do-i-add-time-varying-parameters). +`lateral.river.inflow`: positive for inflows, negative for abstraction. If parameter is time +varying, add it to the correct section, see [below](#how-do-i-add-time-varying-parameters). ### How do I add time-varying parameters? @@ -20,4 +23,6 @@ Either through cyclic (add parameter to `cyclic` list in the toml), or to the fo ### How do I add different output? -See [here for csv output](./toml_file.qmd#output-csv-section), [here for scalar netcdf data](./toml_file.qmd#scalar-data), and [here for gridded output](./toml_file.qmd#output-netcdf-section). \ No newline at end of file +See [here for csv output](./toml_file.qmd#output-csv-section), [here for scalar netcdf +data](./toml_file.qmd#scalar-data), and [here for gridded +output](./toml_file.qmd#output-netcdf-section). \ No newline at end of file diff --git a/docs/user_guide/fews.qmd b/docs/user_guide/fews.qmd index 7a0d095e3..f0c742218 100644 --- a/docs/user_guide/fews.qmd +++ b/docs/user_guide/fews.qmd @@ -3,15 +3,14 @@ title: Run from Delft-FEWS --- Wflow integrates easily as part of an operational system by linking to the -[Delft-FEWS](https://oss.deltares.nl/web/delft-fews/) platform. Delft-FEWS integrates data -and models, and is for example used in many active flood forecasting systems around the -world. +[Delft-FEWS](https://oss.deltares.nl/web/delft-fews/) platform. Delft-FEWS integrates data and +models, and is for example used in many active flood forecasting systems around the world. -This can be done without a model adapter that provides the interface between Delft-FEWS and -an external model (or module). This is possible because time information in the TOML -configuration file is optional and Delft-FEWS can import and export netCDF files. When time -information is left out from the TOML configuration file, the `starttime`, `endtime` and -`timestepsecs` (timestep) of the run is extracted from the netCDF forcing file by wflow. +This can be done without a model adapter that provides the interface between Delft-FEWS and an +external model (or module). This is possible because time information in the TOML configuration +file is optional and Delft-FEWS can import and export netCDF files. When time information is +left out from the TOML configuration file, the `starttime`, `endtime` and `timestepsecs` +(timestep) of the run is extracted from the netCDF forcing file by wflow. To indicate that a wflow model runs from Delft-FEWS, the following setting needs to be specified in the main section of the TOML configuration file: @@ -22,8 +21,8 @@ fews_run = true # optional, default value is false This ensures that wflow offsets the time handling, to meet the expectations of Delft-FEWS. -It also uses a different format for the log file such that each log message takes up only -one line. That meets the [General Adapter +It also uses a different format for the log file such that each log message takes up only one +line. That meets the [General Adapter logFile](https://publicwiki.deltares.nl/display/FEWSDOC/05+General+Adapter+Module#id-05GeneralAdapterModule-logFile) expectations, which then can get parsed with these Delft-FEWS log parsing settings: diff --git a/docs/user_guide/multi_threading.qmd b/docs/user_guide/multi_threading.qmd index ddadf8349..6437ce6fa 100644 --- a/docs/user_guide/multi_threading.qmd +++ b/docs/user_guide/multi_threading.qmd @@ -4,26 +4,25 @@ title: Multithreading ## Using wflow in Julia -Wflow supports multi-threading execution of the wflow\_sbm model that uses the kinematic -wave approach for river, overland and lateral subsurface flow. Both the vertical SBM concept -and the kinematic wave components of this model can run on multiple threads. The optional -[local inertial model for river flow](@ref config_sbm_gwf_lie_river_land) and the optional -[local inertial model for river (1D) and land (2D)](@ref config_sbm_gwf_lie_river_land), -both part of wflow\_sbm, can also run on multiple threads. The threading functionality for -the kinematic wave may also be useful for the wflow\_sbm model [SBM + Groundwater flow](@ref -config_sbm_gwf). The multi-threading functionality in wflow is considered experimental, see -also the following [issue](https://github.com/Deltares/Wflow.jl/issues/139), where an error -was not thrown running code multi-threaded. Because of this we advise to start with running -a wflow model single-threaded (for example during the testing phase of setting up an new -wflow model). +Wflow supports multi-threading execution of the wflow\_sbm model that uses the kinematic wave +approach for river, overland and lateral subsurface flow. Both the vertical SBM concept and the +kinematic wave components of this model can run on multiple threads. The optional [local +inertial model for river flow](@ref config_sbm_gwf_lie_river_land) and the optional [local +inertial model for river (1D) and land (2D)](@ref config_sbm_gwf_lie_river_land), both part of +wflow\_sbm, can also run on multiple threads. The threading functionality for the kinematic +wave may also be useful for the wflow\_sbm model [SBM + Groundwater flow](@ref config_sbm_gwf). +The multi-threading functionality in wflow is considered experimental, see also the following +[issue](https://github.com/Deltares/Wflow.jl/issues/139), where an error was not thrown running +code multi-threaded. Because of this we advise to start with running a wflow model +single-threaded (for example during the testing phase of setting up an new wflow model). For information on how to start Julia with multiple threads we refer to [How to start Julia with multiple threads](https://docs.julialang.org/en/v1/manual/multi-threading/#Starting-Julia-with-multiple-threads). -Additionally, when running Julia + wflow via the command line (note that this is different -from the `wflow_cli`), it is possible to define the number of threads via the `-t` flag. -An example where we start Julia with three threads: +Additionally, when running Julia + wflow via the command line (note that this is different from +the `wflow_cli`), it is possible to define the number of threads via the `-t` flag. An example +where we start Julia with three threads: ```bash julia -t 3 -e 'using Wflow; Wflow.run()' path/to/config.toml @@ -31,11 +30,11 @@ julia -t 3 -e 'using Wflow; Wflow.run()' path/to/config.toml ## Using the command line interface -As explained above, we need to start julia with multiple threads to make use of this -speedup. For `wflow_cli`, the only way to do this is by setting the `JULIA_NUM_THREADS` -environment variable, as explained in [these julia +As explained above, we need to start julia with multiple threads to make use of this speedup. +For `wflow_cli`, the only way to do this is by setting the `JULIA_NUM_THREADS` environment +variable, as explained in [these julia docs](https://docs.julialang.org/en/v1/manual/multi-threading/#Starting-Julia-with-multiple-threads). When a model run starts, among the run information the number of threads that are used is -printed, so `nthreads() = 4` means 4 threads are used, because `JULIA_NUM_THREADS` has been -set to 4. \ No newline at end of file +printed, so `nthreads() = 4` means 4 threads are used, because `JULIA_NUM_THREADS` has been set +to 4. \ No newline at end of file diff --git a/docs/user_guide/output.qmd b/docs/user_guide/output.qmd index b71a24bad..0c08b34cb 100644 --- a/docs/user_guide/output.qmd +++ b/docs/user_guide/output.qmd @@ -5,8 +5,8 @@ title: Handling model output After running the model example from the previous step 4, the model results can be found in `data/output_moselle_simple.csv`. -If required, it is also possible to output netCDF files as output, by modifying the TOML -file. An example is shown below: +If required, it is also possible to output netCDF files as output, by modifying the TOML file. +An example is shown below: ```toml # Spatial output diff --git a/docs/user_guide/required_files.qmd b/docs/user_guide/required_files.qmd index fe6c61941..05cf3986e 100644 --- a/docs/user_guide/required_files.qmd +++ b/docs/user_guide/required_files.qmd @@ -2,24 +2,27 @@ title: Required files --- -In order to run wflow, several files are required. These consist of a settings file and -input data. The input data is typically separated into static maps and forcing data, both -are supplied in a netCDF file, except for lake storage and rating curves that are supplied -via CSV files. A brief overview of the different files: - - - The `settings.toml` file contains information on the simulation period, links to the - input files (and their names in the netCDF files), and links the correct names of the - variables in the netCDF files to the variables and parameters of wflow. - - The `staticmaps.nc` file contains spatial information on the elevation, locations of the - gauges, land-use, drainage direction, etc. This file can also contain maps with parameter +To run wflow, several files are required. These include a settings file and input data. The +input data is typically separated into static maps and forcing data, and both are provided in +netCDF files, except for lake storage and rating curves that are supplied via CSV files. Below +is a brief overview of the different files: + + - The `settings.toml` file contains information on the simulation period, links to the input + files (and their names in the netCDF files), and connect the correct variable names in the + netCDF files to the variables and parameters of wflow. + - The `staticmaps.nc` file contains spatial information such as elevation, gauge locations, + land use, and drainage direction, etc. This file can also contain maps with parameter values. - - The `forcing.nc` file contains the precipitation, temperature and potential evaporation - time series (as a 3D array). + - The `forcing.nc` file contains time series data for precipitation, temperature and + potential evaporation (as a 3D array). ## The configuration file (`settings.toml`) -The configuration file contains all relevant settings for running wflow, such as the simulation period, the model settings, the mapping between input files and (internal) model parameters. More details and explanations can be found [here](./toml_file.qmd). An example configuration file is presented below. +The configuration file contains all relevant settings for running wflow, such as the simulation +period, the model settings, the mapping between input files and (internal) model parameters. +More details and explanations can be found [here](./toml_file.qmd). An example configuration +file is presented below.
Click to show example .toml file @@ -33,8 +36,8 @@ The configuration file contains all relevant settings for running wflow, such as The list below contains a brief overview of several essential static maps required to run wflow. These NC variables names refer to the example data of the wflow\_sbm + kinematic wave -model (see [here](@ref wflow_sbm_data)). Example data for the other model configurations can -be found [here](@ref sample_data). +model (see [here](@ref wflow_sbm_data)). Example data for the other model configurations can be +found [here](@ref sample_data). Description | NC variable name | unit --- | --- | --- @@ -47,9 +50,9 @@ Land slope | `Slope` | m m$^{-1}$ River slope | `RiverSlope` | m m$^{-1}$ As mentioned before, the model parameters can also be defined as spatial maps. They can be -included in the same netCDF file, as long as their variable names are correctly mapped in -the TOML settings file. See the section on [example models](@ref sample_data) on how to -use this functionality. +included in the same netCDF file, as long as their variable names are correctly mapped in the +TOML settings file. See the section on [example models](@ref sample_data) on how to use this +functionality. ::: callout-important When using cyclic data, three different options are supported: @@ -57,7 +60,10 @@ When using cyclic data, three different options are supported: - 365 (daily, where Feb. 29 is not present, so the value for Feb. 28 is taken. Dec. 31 is represented by DOY 365) - 366 (where Feb. 29 represents DOY 60, and Dec. 31 DOY 366) -In contrast to the right-labelling of the forcing data (see below), the DOY/month of the current timestep is used. For example: to simulate 2023-06-14 00:00:00 (with a daily timestep), the DOY value at position 6 (when 12 values are provided), 165 (when 365 values are provided) or 166 (when 366 values are provided). +In contrast to the right-labelling of the forcing data (see below), the DOY/month of the +current time step is used. For example: to simulate 2023-06-14 00:00:00 (with a daily time +step), the DOY value at position 6 (when 12 values are provided), 165 (when 365 values are +provided) or 166 (when 366 values are provided). ::: ## The forcing data (`forcing.nc`) @@ -68,12 +74,12 @@ precipitation, temperature and potential evaporation. The code snippet below sho of the example file (downloaded [here](@ref wflow_sbm_data)), and displaying the content with `NCDatasets` in Julia. As can be seen, each forcing variable (`precip`, `pet` and `temp`) consists of a three-dimensional dataset (`x`, `y`, and `time`), and each timestep consists of a -two-dimensional map with values at each gridcell. Only values within the basin are required. +two-dimensional map with values at each grid cell. Only values within the basin are required. ::: callout-important -Wflow expects the foring to be "right-labelled". This means that e.g. daily precipitation -at 2023-06-15 00:00:00 is the accumulated total precipitation between 2023-06-14 -00:00:00 and 2023-06-15 00:00:00. +Wflow expects the forcing to be "right-labelled". This means that e.g. daily precipitation at +2023-06-15 00:00:00 is the accumulated total precipitation between 2023-06-14 00:00:00 and +2023-06-15 00:00:00. ::: ::: callout-note diff --git a/docs/user_guide/toml_file.qmd b/docs/user_guide/toml_file.qmd index 8a922e582..1c7211d6d 100644 --- a/docs/user_guide/toml_file.qmd +++ b/docs/user_guide/toml_file.qmd @@ -3,26 +3,25 @@ title: Simulation settings --- All model configuration and settings are passed through the .toml file, which contains all -relevant settings, information on the model configuration, simulation period, input files, and -parameters. The settings are provided in a TOML file. The settings file is structured in -several sections, which are explained below. The filepaths that are provided in this file are -relative to the location of the TOML file, or to `dir_input` and `dir_output` if they are -given. +relevant settings, information about the model configuration, simulation period, input files, +and parameters. The settings are provided in a TOML file. The settings file is structured in +several sections, which are explained below. The file paths provided in this file are relative +to the location of the TOML file, or to `dir_input` and `dir_output` if they are specified. ## General time info -Time information is optional. When omitted, wflow will perform computations for each -timestamp in the forcing netCDF file, except for the first forcing timestamp, which is -considered equal to the initial conditions of the wflow model (state time). If you wish to -calculate a subset of this time range, or a different timestep, you can specify a -`starttime`, `endtime` and `timestepsecs`. The `starttime` is defined as the model state -time. In the TOML file settings below, the `starttime` is 2000-01-01T00:00:00 (state time) -and the first update (and output) of the wflow model is at 2000-01-02T00:00:00. The -`time_units` optional information is used by the `writer` of the model, for model output in -netCDF format. The `calendar` option allows you to calculate in one of the different [CF -conventions calendars](http://cfconventions.org/cf-conventions/cf-conventions.html#calendar) -provided by the [CFTime.jl package](https://juliageo.org/CFTime.jl/latest/), such as -`"360_day"`. This is useful if you want to calculate climate scenarios which are sometimes -provided in these alternative calendars. +Time information is optional. When omitted, wflow will perform computations for each timestamp +in the forcing netCDF file, except for the first forcing timestamp, which is considered equal +to the initial conditions of the wflow model (state time). If you wish to calculate a subset of +this time range, or a different timestep, you can specify a `starttime`, `endtime` and +`timestepsecs`. The `starttime` is defined as the model state time. In the TOML file settings +below, the `starttime` is 2000-01-01T00:00:00 (state time) and the first update (and output) of +the wflow model is at 2000-01-02T00:00:00. The `time_units` optional information is used by the +`writer` of the model, for model output in netCDF format. The `calendar` option allows you to +calculate in one of the different [CF conventions +calendars](http://cfconventions.org/cf-conventions/cf-conventions.html#calendar) provided by +the [CFTime.jl package](https://juliageo.org/CFTime.jl/latest/), such as `"360_day"`. This is +useful if you want to calculate climate scenarios which are sometimes provided in these +alternative calendars. ```toml calendar = "standard" # optional, this is the default value @@ -35,11 +34,11 @@ dir_output = "data/output" # optional, default is the path ``` ## Logging -Wflow emits logging messages at various levels such as debug, info, and error. These get -sent to both the terminal as well as a log file. Note that logging to a file is only part of -the `Wflow.run(tomlpath::AbstractString)` method. If you want to debug an issue it can be -helpful to set `loglevel = "debug"` in the TOML. To avoid flooding the screen, debug -messages are only sent to the log file. The following settings will affect the logging: +Wflow prints logging messages at various levels such as debug, info, and error. These messages +are sent to both the terminal and a log file. Note that logging to a file is only part of the +`Wflow.run(tomlpath::AbstractString)` method. If you want to debug an issue, it can be helpful +to set `loglevel = "debug"` in the TOML. To avoid flooding the screen, debug messages are only +sent to the log file. The following settings will affect the logging: ```toml silent = false # optional, default is "false" @@ -48,13 +47,13 @@ path_log = "log.txt" # optional, default is "log.txt" fews_run = false # optional, default value is false ``` -`silent` avoids logging to the terminal, and only writes to the log file. `loglevel` -controls which levels are filtered out; for instance, the default setting `"info"` does not -print any debug-level messages. Note that for finer control, you can also pass an integer -log level. For details, see Julia's -[Logging](https://docs.julialang.org/en/v1/stdlib/Logging/#Log-event-structure) -documentation. `path_log` sets the desired output path for the log file. For information on -`fews_run`, see [Run from Delft-FEWS](@ref run_fews). +`silent` avoids logging to the terminal, and only writes to the log file. `loglevel` controls +which levels are filtered out; for instance, the default setting `"info"` does not print any +debug-level messages. Note that for finer control, you can also pass an integer log level. For +details, see Julia's +[Logging](https://docs.julialang.org/en/v1/stdlib/Logging/#Log-event-structure) documentation. +`path_log` sets the desired output path for the log file. For information on `fews_run`, see +[Run from Delft-FEWS](@ref run_fews). ## Model section Model-specific settings can be included in the model section of the TOML file. @@ -74,12 +73,12 @@ min_streamorder_land = 4 # minimum stream order to delineate subbasin ## State options The `state` section in the TOML file provides information about the location of input and -output states of the model. This section is mostly relevant if the model needs to start with -a "warm" state (i.e. based on the results of a previous simulation). The example below shows -how to save the output states of the current simulation, so it can be used to initialize -another model in the future. Details on the settings required to start a model with a warm -state can be found in the [additional model options](@ref reinit). If it is not required to -store the outstates of the current simulation, the entire `state` section can be removed. +output states of the model. This section is mostly relevant if the model needs to start with a +"warm" state (i.e. based on the results of a previous simulation). The example below shows how +to save the output states of the current simulation, so it can be used to initialize another +model in the future. Details on the settings required to start a model with a warm state can be +found in the [additional model options](@ref reinit). If it is not required to store the +outstates of the current simulation, the entire `state` section can be removed. ```toml [state] @@ -113,16 +112,15 @@ h_av = "h_av_land" ## Input section The `input` section of the TOML file contains information about the input forcing and model -parameters files (in netCDF format). Forcing is applied to the vertical component of the -model, and needs to be mapped to the external netCDF variable name. `forcing` lists the -internal model forcing parameters, and these are mapped to the external netCDF variables -listed under the section `[input.vertical]`. It is possible to provide cyclic parameters to -the model (minimum time step of 1 day). In the example below, the internal -`vertical.leaf_area_index` model parameter is mapped to the external netCDF variable "LAI" -variable. Cyclic time inputs of parameters can be different (e.g., daily or monthly). The -`time` dimension name of these cylic input parameters in the model parameter netCDF file -should start with "time". If a model parameter is not mapped, a default value will be used, -if available. +parameters files (in netCDF format). Forcing is applied to the vertical component of the model, +and needs to be mapped to the external netCDF variable name. `forcing` lists the internal model +forcing parameters, and these are mapped to the external netCDF variables listed under the +section `[input.vertical]`. It is possible to provide cyclic parameters to the model (minimum +time step of 1 day). In the example below, the internal `vertical.leaf_area_index` model +parameter is mapped to the external netCDF variable "LAI" variable. Cyclic time inputs of +parameters can be different (e.g., daily or monthly). The `time` dimension name of these cylic +input parameters in the model parameter netCDF file should start with "time". If a model +parameter is not mapped, a default value will be used, if available. ```toml [input] @@ -236,23 +234,25 @@ h = "h_land" ``` ### Scalar data -Besides gridded data, it is also possible to write scalar data to a netCDF file. Below is an -example that writes scalar data to the file "output\_scalar\_moselle.nc". For each netCDF -variable a `name` (external variable name) and `parameter` (internal model parameter) is -required. A `reducer` can be specified to apply to the model output, see for more -information the following section [Output CSV section](@ref). When a `map` is provided to -extract data for certain locations (e.g. `gauges`) or areas (e.g. `subcatchment`), the -netCDF location names are extracted from these maps. For a specific location (grid cell) a -`location` is required. For layered model parameters and variables that have an extra -dimension `layer` and are part of the vertical `sbm` concept it is possible to specify an -internal layer index (see also example below). If multiple layers are desired, this can be -specified in separate `[[netcdf.variable]]` entries. Note that the specification of the -extra dimension is not optional when wflow is integrated with Delft-FEWS, for netCDF scalar -data an extra dimension is not allowed by the `importNetcdfActivity` of the Delft-FEWS -General Adapter. In the section [Output CSV section](@ref), similar functionality is -available for CSV. For integration with Delft-FEWS, see also [Run from Delft-FEWS](@ref -run_fews), it is recommended to write scalar data to netCDF format since the General Adapter -of Delft-FEWS can ingest this data format directly. +In addition to gridded data, scalar data can also be written to a netCDF file. Below is an +example that shows how to write scalar data to the file "output\_scalar\_moselle.nc". For each +netCDF variable, a `name` (external variable name) and a `parameter` (internal model parameter) +are required. A `reducer` can be specified to apply to the model output. See more details in +the [Output CSV section](@ref) section. If a `map` is provided to extract data for specific +locations (e.g. `gauges`) or areas (e.g. `subcatchment`), the netCDF location names are +extracted from these maps. For a specific location (grid cell) a `location` is required. For +layered model parameters and variables that have an extra `layer` dimension and are part of the +vertical `sbm` concept, an internal layer index can be specified (an example is provided +below). Similarly, for model parameters and variables that have an extra dimension `classes` +and are part of the vertical `FLEXTopo` concept, the class name can be specified. If multiple +layers or classes are desired, this can be specified in separate `[[netcdf.variable]]` entries. +Note that the additional dimension should be specified when wflow is integrated with +Delft-FEWS, for netCDF scalar data an extra dimension is not allowed by the +`importNetcdfActivity` of the Delft-FEWS General Adapter. In the section [Output CSV +section](@ref), similar functionality is available for CSV. For integration with Delft-FEWS, it +is recommended to write scalar data to netCDF format, as the General Adapter of Delft-FEWS can +directly ingest data from netCDF files. For more information, see [Run from Delft-FEWS](@ref +run_fews). ```toml [netcdf] @@ -280,10 +280,10 @@ parameter = "vertical.temperature" ``` ## Output CSV section -Model output can also be written to a CSV file. Below is an example that writes model output -to the file "output_moselle.csv". For each CSV column, a `header` and `parameter` (internal -model parameter) are required. A `reducer` can be specified to apply to the model output, -with the following available reducers: +Model output can also be written to a CSV file. Below is an example that writes model output to +the file "output_moselle.csv". For each CSV column, a `header` and `parameter` (internal model +parameter) are required. A `reducer` can be specified to apply to the model output, with the +following available reducers: + maximum + minimum @@ -294,15 +294,17 @@ with the following available reducers: + last + only -with `only` as the default. To extract data for a specific location (grid cell), the `index` -of the vector, the coordinates `coordinate.x` and `coordinate.y`, or the x and y indices of -the 2D array (`index.x` and `index.y`) can be provided. Finally a `map` can be provided to +with `only` as the default. To extract data for a specific location (grid cell), the `index` of +the vector, the coordinates `coordinate.x` and `coordinate.y`, or the x and y indices of the 2D +array (`index.x` and `index.y`) can be provided. Additionally, a `map` can be provided to extract data for certain locations (e.g. `gauges`) or areas (e.g. `subcatchment`). In this -case a single entry can lead to multiple columns in the CSV file, which will be of the form +case, a single entry can lead to multiple columns in the CSV file, which will be of the form `header_id`, e.g. `Q_20`, for a gauge with integer ID 20. For layered model parameters and variables that have an extra dimension `layer` and are part of the vertical `sbm` concept an -internal layer index (see also example below) should be specified. If multiple layers are -desired, this can be specified in separate `[[csv.column]]` entries. +internal layer index (see also example below) should be specified. For model parameters and +variables that have an extra dimension `classes` and are part of the vertical `FLEXTopo` +concept, it is possible to specify the class name. If multiple layers or classes are desired, +this can be specified in separate `[[csv.column]]` entries. The double brackets in `[[csv.column]]` follow TOML syntax, indicating that it is part of a list. You can specify as many entries as you want. @@ -369,8 +371,8 @@ waterfrac = "WaterFrac" cfmax.value = 2.5 ``` -For input parameters with an extra dimension (e.g. `layer` or `classes`), one uniform value -can be provided or a list of values that matches the length of the additional dimension. For +For input parameters with an extra dimension (e.g. `layer` or `classes`), one uniform value can +be provided or a list of values that matches the length of the additional dimension. For example, a list of values can be provided for input parameter `c` as follows: ```toml @@ -380,8 +382,8 @@ waterfrac = "WaterFrac" c.value = [10.5, 11.25, 9.5, 7.0] ``` -To change the forcing variable `precipitation` with a `scale` factor of 1.5 and -an `offset` of 0.5: +To change the forcing variable `precipitation` with a `scale` factor of 1.5 and an `offset` of +0.5: ```toml [input.vertical.precipitation] @@ -391,9 +393,9 @@ offset = 0.5 ``` For input parameters with an extra dimension, it is also possible to modify multiple indices -simultaneously with different `scale` and `offset` values. In the example below, the -external netCDF variable `c` is modified at `layer` index 1 and 2, with a `scale` factor of -2.0 and 1.5 respectively, and an `offset` of 0.0 for both indices: +simultaneously with different `scale` and `offset` values. In the example below, the external +netCDF variable `c` is modified at `layer` index 1 and 2, with a `scale` factor of 2.0 and 1.5 +respectively, and an `offset` of 0.0 for both indices: ```toml [input.vertical.c] @@ -404,9 +406,9 @@ layer = [1, 2] ``` ## Fixed forcing values -It is possible to set fixed values for forcing parameters through the TOML file. For -example, to set `temperature` to a fixed value of 10 ``\degree``C, the complete `forcing` -list is required: +It is possible to set fixed values for forcing parameters through the TOML file. For example, +to set `temperature` to a fixed value of 10 ``\degree``C, the complete `forcing` list is +required: ```toml forcing = [ diff --git a/docs/user_guide/warm_states.qmd b/docs/user_guide/warm_states.qmd index 35faac9b6..28cdd2170 100644 --- a/docs/user_guide/warm_states.qmd +++ b/docs/user_guide/warm_states.qmd @@ -3,18 +3,17 @@ title: Starting with "warm" states --- The `state` section in the TOML file provides information on the input file if the model is -initialized with a warm state (`path_input`) and to what file the states are written at the -end of the model run (`path_output`). Please note that the model setting `reinit` needs to -be set to `false` in order to initialize the model with states from the file located at -`path_input`. A mapping between external state names and internal model states is required. -This information is specified for each model component, the `vertical` model and `lateral` -model components. In the example below the `vertical` component represents the SBM concept, -and for the `lateral` components there is a `river` (including optional `reservoir`, `lake` -and `floodplain` components), `land` and `subsurface` domain. The internal model states are -listed on the left side, and the external state names are listed on the right side. Note -that `path_input` is only required when `reinit` is set to false. `path_output` is optional, -an output state file is only written when it is defined. If neither is set, the entire -`state` section can be left out. +initialized with a warm state (`path_input`) and to what file the states are written at the end +of the model run (`path_output`). Please note that the model setting `reinit` needs to be set +to `false` in order to initialize the model with states from the file located at `path_input`. +A mapping between external state names and internal model states is required. This information +is specified for each model component, the `vertical` model and `lateral` model components. In +the example below the `vertical` component represents the SBM concept, and for the `lateral` +components there is a `river` (including optional `reservoir`, `lake` and `floodplain` +components), `land` and `subsurface` domain. The internal model states are listed on the left +side, and the external state names are listed on the right side. Note that `path_input` is only +required when `reinit` is set to false. `path_output` is optional, an output state file is only +written when it is defined. If neither is set, the entire `state` section can be left out. ```toml [model] diff --git a/docs/user_guide/zmq_server.qmd b/docs/user_guide/zmq_server.qmd index 453279755..06f6595ed 100644 --- a/docs/user_guide/zmq_server.qmd +++ b/docs/user_guide/zmq_server.qmd @@ -3,7 +3,7 @@ title: Run wflow as a ZMQ Server --- It is possible to run wflow as a ZMQ Server, for example for the coupling to the -[OpenDA](https://openda.org/) software for data-assimilation. The code for the wflow ZMQ -Server is not part of the Wflow.jl package, and is located +[OpenDA](https://openda.org/) software for data-assimilation. The code for the wflow ZMQ Server +is not part of the Wflow.jl package, and is located [here](https://github.com/Deltares/Wflow.jl/tree/master/server).