From 575eca97c577f92d29d938f2f7b60373ece81884 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Sat, 26 Oct 2024 19:55:17 +0200 Subject: [PATCH] Update formulas and units in docs for lateral concepts and sbm --- .gitignore | 1 + docs/make.jl | 7 + docs/src/model_docs/lateral/gwf.md | 85 ++++--- docs/src/model_docs/lateral/kinwave.md | 65 ++--- docs/src/model_docs/lateral/local-inertial.md | 61 ++--- docs/src/model_docs/lateral/sediment_flux.md | 231 ++++++++--------- docs/src/model_docs/lateral/waterbodies.md | 54 ++-- docs/src/model_docs/vertical/sbm.md | 234 +++++++++--------- 8 files changed, 378 insertions(+), 360 deletions(-) diff --git a/.gitignore b/.gitignore index 604242d59..9d64d76fc 100644 --- a/.gitignore +++ b/.gitignore @@ -13,6 +13,7 @@ # Document generation /docs/build/ /docs/site/ +/docs/Manifest.toml # Temporary files /dev/ diff --git a/docs/make.jl b/docs/make.jl index a710b5fb3..25d1d749e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -53,6 +53,13 @@ makedocs(; canonical = "https://deltares.github.io/Wflow.jl", assets = String[], collapselevel = 2, + mathengine = Documenter.KaTeX(Dict( + :macros => Dict( + raw"\SI" => raw"{#1\;\mathrm{#2}}", + raw"\SIb" => raw"{#1\;[\mathrm{#2}}]", # b for brackets + raw"\subtext" => raw"#1_\text{#2}", + ) + )) ), pages = pages, ) diff --git a/docs/src/model_docs/lateral/gwf.md b/docs/src/model_docs/lateral/gwf.md index 8f932c8a0..5914ed81e 100644 --- a/docs/src/model_docs/lateral/gwf.md +++ b/docs/src/model_docs/lateral/gwf.md @@ -11,40 +11,39 @@ Single layer groundwater flow requires the four following components, and each i 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 ``H`` [m] over the aquifer and transmissivity -``kH`` [m``^2`` d``^{-1}``] is a constant (``k`` [m d``^{-1}``] is the horizontal hydraulic +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 1e-5 (stiff) and 0.01 +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). +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: ```math - S \frac{\phi}{\delta t} = \frac{\delta}{\delta x} (kH \frac{\phi}{\delta x}) + Q + S \frac{\partial \phi}{\partial t} = \frac{\partial}{\partial x} \left(kH \frac{\phi}{\delta x}\right) + Q ``` -where ``S`` [m m``^{-1}``] is storativity (or specific yield), ``\phi`` [m] is hydraulic -head, ``t`` is time, ``k`` [m t``^{-1}``] is horizontal hydraulic conductivity, ``H`` [m] is -the (saturated) aquifer height: groundwater level - aquifer bottom elevation and ``Q`` [m -t``^{-1}``] represents fluxes from boundary conditions (e.g. recharge or abstraction), see +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 simplest finite difference formulation is forward in time, central in space, and can be written as: ```math - 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_ᵢ + 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}`` +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. @@ -55,17 +54,17 @@ Conductance ``C`` is defined as: C = \frac{kH w}{l} ``` -where ``w`` [m] is the width of the cell to cell connection, and ``l`` [m] is the length of +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: ```math - 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)} + 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 ``H`` [m] is the aquifer top - aquifer bottom, and ``k``, ``l_i`` is the length in +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 @@ -89,12 +88,12 @@ Finally, a stable time step size can be computed given the forward-in-time, cent scheme, based on the following criterion from Chu and Willis (1984): ```math -\frac{\Delta t k H}{(\Delta x \Delta y S)} \le \frac{1}{4} +\frac{\Delta t k H}{\Delta x \Delta y S} \le \frac{1}{4} ``` -where ``\Delta t`` [d] is the stable time step size, ``\Delta x`` [m] is the cell length in -the ``x`` direction and ``\Delta y`` [m] is the cell length in the ``y`` direction, ``k`` is -the horizontal hydraulic conductivity [m``^2`` d``^{-1}``] and ``H`` [m] is the saturated -thickness of the aquifer. For each cell ``\frac{(\Delta x \Delta y S)}{k H}`` is +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. @@ -146,38 +145,44 @@ The flux between river and aquifer is calculated using Darcy's law following the MODFLOW: ```math - Q_{riv} = \Bigg\lbrace{C_{i} \,\text{min}(h_{riv} - B_{riv}, h_{riv} - \phi), \,h_{riv} > \phi \atop C_{e} (h_{riv} - \phi) , \,h_{riv} \leq \phi} + \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 \\ + C_e (\subtext{h}{riv} - \phi) &\text{ if }\quad \subtext{h}{riv} \le \phi + \end{cases} + \end{align*} ``` -where ``Q_{riv}`` is the exchange flux from river to aquifer [L``^3`` T``^{-1}``], ``C_i`` -[L``^2`` T``^{-1}``] is the river bed infiltration conductance, ``C_e`` [L``^2`` T``^{-1}``] -is the river bed exfiltration conductance, ``B_{riv}`` the bottom of the river bed [L], -``h_{riv}`` is the river stage [L] and ``\phi`` is the hydraulic head in the river cell [L]. + +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 exchange flux (river to aquifer) ``Q_{riv}`` is an output variable (field `flux` of the +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 `h` [m] of the river kinematic wave in combination with +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. ### Drainage The flux from drains to the aquifer is calculated as follows: ```math -Q_{drain} = C_{drain} \text{min}(0, h_{drain} - \phi) +\subtext{Q}{drain} = \subtext{C}{drain} \min(0, \subtext{h}{drain} - \phi) ``` -where ``Q_{drain}`` is the exchange flux from drains to aquifer [L``^3`` T``^{-1}``], -``C_{drain}`` [L``^2`` T``^{-1}``] is the drain conductance, ``h_{drain}`` is the drain -elevation [L] and ``\phi`` is the hydraulic head in the cell with drainage [L]. +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. -The Table in the Groundwater flow [drainage boundary condition](@ref gwf_drainage_params) +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) ``Q_{drain}`` is an output variable (field `flux` of +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 @@ -194,10 +199,10 @@ The recharge flux ``Q_{r}`` to the aquifer is calculated as follows: ```math Q_{r} = R \, A ``` -with ``R`` the recharge rate [L T``^{-1}``] and ``A`` the area [L``^2`` ] of the aquifer +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) +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. @@ -212,15 +217,15 @@ multiplied by the `area` field of the aquifer. 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 ``Q_{hb}`` [L``^3`` T``^{-1}``] is calculated as follows: +the boundary ``\SIb{Q_{hb}}{L^3 T^{-1}}`` is calculated as follows: ```math Q_{hb} = C_{hb} (\phi_{hb} - \phi) ``` -with ``C_{hb}`` the conductance of the head boundary [L``^2`` T``^{-1}``], ``\phi_{hb}`` the -head [L] of the head boundary and ``\phi`` the head of the aquifer cell. +with ``\SIb{C_{hb}}{L^2 T^{-1}}`` the conductance of the head boundary, ``\SIb{\phi_{hb}}{L}`` the +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) +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 @@ -232,12 +237,12 @@ value. This boundary is not (yet) part of the model `SBM + Groundwater flow`. ### Well boundary -A volumetric well rate [L``^3`` T``^{-1}``] can be specified as a boundary condition. +A volumetric well rate ``\SIb{}{L^3 T^{-1}}`` can be specified as a boundary condition. 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 ``Q_{well}`` can be can be specified as a fixed or time dependent +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). diff --git a/docs/src/model_docs/lateral/kinwave.md b/docs/src/model_docs/lateral/kinwave.md index d0e744e73..33fa63705 100644 --- a/docs/src/model_docs/lateral/kinwave.md +++ b/docs/src/model_docs/lateral/kinwave.md @@ -1,23 +1,24 @@ # [Kinematic wave] (@id kin_wave) ## Surface routing -The main flow routing scheme available in Wflow.jl is the kinematic wave approach for -channel and overland flow, assuming that the topography controls water flow mostly. The +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): ```math - \dfrac{dQ}{dx} + \dfrac{dA}{dt} = q \\~\\ - A = \alpha Q^{\beta} + \dfrac{\partial Q}{\partial x} + \dfrac{\partial A}{\partial t} = q, \\~\\ + A = \alpha Q^{\beta}. ``` These equations can then be combined as a function of streamflow only: ```math - \dfrac{dQ}{dx} + \alpha \beta Q^{\beta - 1} \dfrac{dQ}{dt} = q + \dfrac{\partial Q}{\partial x} + \alpha \beta Q^{\beta - 1} \dfrac{\partial Q}{\partial t} = q. ``` -where ``Q`` is the surface runoff in the kinematic wave [m``^3``/s], ``x`` is the length of -the runoff pathway [m], ``A`` is the cross-section area of the runoff pathway [m``^{2}``], -``t`` is the integration timestep [s] and ``\alpha`` and ``\beta`` are coefficients. + +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 +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) @@ -36,7 +37,7 @@ 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 1.0 m) for the river can be provided as follows: +is ``\SI{1.0}{m}``) for the river can be provided as follows: ```toml [input.lateral.river] @@ -50,12 +51,12 @@ Simplified [reservoir and lake](@ref reservoir_lake) models can be included as p river kinematic wave network. ## Inflow -External water (supply/abstraction) `inflow` [m``^3`` s``^{-1}``] can be added to the +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)). ## Abstractions -Abstractions from the river through the variable `abstraction` [m``^3`` s``{-1}``] are +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 @@ -72,29 +73,30 @@ detailed description) are available: - `layered` - `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 ``w`` [m] +drained laterally by saturated downslope subsurface flow for a slope with width ``\SIb{w}{m}`` according to: ```math - Q = \begin{cases} - \frac{K_0\tan(\beta)}{f}\left(e^{(-fz_{i})}-e^{(-fz_\mathrm{exp})}\right) w + - K_0e^{(-fz_\mathrm{exp})}(z_t-z_\mathrm{exp})\tan(\beta) w & \text{if $z_i < z_\mathrm{exp}$}\\ + Q = K_0\tan(\beta)w\begin{cases} + \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}$}\\ \\ - K_0e^{(-fz_\mathrm{exp})}(z_t - z_i)\tan(\beta) w & \text{if $z_i \ge z_\mathrm{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, ``Q`` is subsurface flow [m``^{3}`` d``^{-1}``], -``K_0`` is the saturated hydraulic conductivity at the soil surface [m d``^{-1}``], ``z_i`` -is the water table depth [m], ``z_{t}`` is the total soil depth [m], ``f`` is a scaling -parameter [m``^{-1}``] that controls the decrease of ``K_0`` with depth and -``z_\mathrm{exp}`` [m] is the depth from soil surface for which the exponential decline of -``K_0`` is valid. For the `exponential` profile, ``z_\mathrm{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: ```math (\theta_s-\theta_r)w\frac{\partial h}{\partial t} = -\frac{\partial Q}{\partial x} + wr ``` -where ``h`` is the water table height [m], ``x`` is the distance downslope [m], and ``r`` is -the net input rate [m d``^{-1}``] 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: ```math \frac{\partial Q}{\partial t} = -c\frac{\partial Q}{\partial x} + cwr @@ -102,23 +104,22 @@ the net input rate [m d``^{-1}``] to the saturated store. Substituting for ``h where celerity ``c`` is calculated as follows: ```math - c = \begin{cases} - \frac{K_0e^{(-fz_{i})}\tan(\beta)}{(\theta_s-\theta_r)} - + \frac{K_0e^{(-fz_\mathrm{exp})}\tan(\beta)}{(\theta_s-\theta_r)} & \text{if $z_i < z_\mathrm{exp}$}\\ + c = \frac{K_0 \tan(\beta)}{\theta_s-\theta_r}\begin{cases} + e^{-fz_i} + + e^{-f\subtext{z}{exp}} & \text{if $z_i < \subtext{z}{exp}$}\\ \\ - \frac{K_0e^{(-fz_\mathrm{exp})}\tan(\beta)}{(\theta_s-\theta_r)} & \text{if $z_i \ge z_\mathrm{exp}$}. + e^{-f\subtext{z}{exp}} & \text{if $z_i \ge \subtext{z}{exp}$}. \end{cases} ``` For the `layered` and `layered_exponential` profiles the equivalent horizontal hydraulic -conductivity ``K_h`` [m d``^{-1}``] is calculated for water table height ``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: ```math Q = K_h h \tan(\beta) w, ``` and celerity ``c`` is given by: ```math - c = \frac{K_h \tan(\beta)}{(\theta_s-\theta_r)}. + c = \frac{K_h \tan(\beta)}{\theta_s-\theta_r}. ``` The kinematic wave equation for lateral subsurface flow is solved iteratively using Newton's diff --git a/docs/src/model_docs/lateral/local-inertial.md b/docs/src/model_docs/lateral/local-inertial.md index ed320c754..75c881ee1 100644 --- a/docs/src/model_docs/lateral/local-inertial.md +++ b/docs/src/model_docs/lateral/local-inertial.md @@ -8,33 +8,33 @@ of the local inertial approximation on a staggered grid is as follows (Bates et ```math 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 ``Q_{t+\Delta t}`` is the river flow [m``^3``/s] at time step ``t+\Delta t``, ``g`` is -acceleration due to gravity [m/s``^2``], ``A_t`` is the cross sectional flow area at the -previous time step, ``R_t`` is the hydraulic radius at the previous time step, ``Q_t`` is -the river flow [m``^3``/s] at the previous time step, ``S_t`` is the water surface slope at -the previous time step and ``n`` is the Manning's roughness coefficient [m``^{-1/3}`` s]. +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: ```math -h^{t+\Delta t} = h^t + \Delta t \frac{Q^{t+\Delta t}_{src} - Q^{t+\Delta t}_{dst}}{A} +h^{t+\Delta t} = h^t + \Delta t \frac{\subtext{Q^{t+\Delta t}}{src} - \subtext{Q^{t+\Delta t}}{dst}}{A} ``` -where ``h^{t+\Delta t}`` is the water depth [m] at time step ``t+\Delta t``, ``h^t`` is the -water depth [m] at the previous time step, ``A`` is the river area [m``^2``] and ``Q_{src}`` -and ``Q_{dst}`` represent river flow [m``^3``/s] at the upstream and downstream link of the +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): ```math -\Delta t = min(\alpha \frac{\Delta x_i}{\sqrt{(gh_i)}}) +\Delta t = \alpha \min_i\left(\frac{\Delta x_i}{\sqrt{gh_i}}\right) ``` -where ``\sqrt{(gh_i)}`` is the wave celerity for river cell ``i`` , ``\Delta x_i`` is the -river length [m] for river cell ``i`` and ``\alpha`` is a coefficient (typically between 0.2 -and 0.7) to enhance the stability of the simulation. +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. In the TOML file the following properties related to the local inertial model can be provided for the `sbm` and `sbm_gwf` model types: @@ -67,24 +67,24 @@ The momentum equation is most stable for low slope environments, and to keep the 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 ``Fr`` on a link +et al. (2017) in the Landlab v1.0 OverlandFlow component. The froude number ``\mathrm{Fr}`` on a link is calculated as follows: ```math - Fr = \frac{u}{\sqrt{(gh_f)}} + \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 +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. +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 width, river bed elevation and Manning's roughness coefficient are copied from the upstream -river cell. The river length [m] of the boundary cell can be set through the TOML file with -`riverlength_bc`, and has a default value of 10 km. The water depth at the boundary cell is -fixed at 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 @@ -106,14 +106,15 @@ Q_{i-1/2}^{n+1} = \frac{\left[ \theta Q_{i-1/2}^{n} +\frac{(1-\theta)}{2}(Q_{(i- n^2 |Q_{i-1/2}^{n}|/(h_f^{7/3} \Delta y)} ``` + where subscripts ``i`` and ``n`` refer to space and time indices, respectively. Subscript -``i-1/2`` is to the link between node ``i`` and ``i-1``, subscript ``i+1/2`` is the link -between node ``i`` and node ``i+1``, and subscript ``i-3/2`` is the link between node ``i-1`` -and node ``i-2``. ``Q`` is the water discharge [m``^3`` s``^{-1}``], ``\eta`` is the water -surface elevation [m], ``h_f`` [m] is the water depth between cells, ``n`` is the Manning's -roughness coefficient [m``^{-1/3}`` s], ``g`` is acceleration due to gravity [m/s``^2``], -``\Delta t`` [s] is the adaptive model time step, ``\Delta x`` [m] is the distance between -two cells and ``\Delta y`` [m] is the flow width. Below the staggered grid and variables of +``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): ![numerical_scheme_almeida](../../images/numerical_scheme_almeida.png) @@ -142,12 +143,12 @@ routing as well as 2D overland flow. The properties `inertial_flow_alpha` and in the [River and floodplain routing](@ref) section of the local inertial model. ## Inflow -External water (supply/abstraction) `inflow` [m``^3`` s``^{-1}``] can be added to the local +External water (supply/abstraction) `inflow` ``\SIb{}{m^3 s^{-1}}`` can be added to the local inertial model for river flow (1D) and river and overland flow combined (1D-2D), as a cyclic parameter or as part of forcing (see also [Input section](@ref)). ## Abstractions -Abstractions from the river through the variable `abstraction` [m``^3`` s``{-1}``] are +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. diff --git a/docs/src/model_docs/lateral/sediment_flux.md b/docs/src/model_docs/lateral/sediment_flux.md index 9a123070e..8ab783ce2 100644 --- a/docs/src/model_docs/lateral/sediment_flux.md +++ b/docs/src/model_docs/lateral/sediment_flux.md @@ -19,37 +19,44 @@ issues, the Yalin transport equation was chosen as it can handle particle differ with no particle differentiation). For land cells, wflow\_sediment assumes that erosion can mobilize 5 classes of sediment: -- Clay (mean diameter of 2 ``\mu``m) -- Silt (mean diameter of 10 ``\mu``m) -- Sand (mean diameter of 200 ``\mu``m) -- Small aggregates (mean diameter of 30 ``\mu``m) -- Large aggregates (mean diameter of 500 ``\mu``m). +- Clay (mean diameter of ``\SI{2}{\mu m}``) +- Silt (mean diameter of ``\SI{10}{\mu m}``) +- Sand (mean diameter of ``\SI{200}{\mu m}``) +- Small aggregates (mean diameter of ``\SI{30}{\mu m}``) +- Large aggregates (mean diameter of ``\SI{50}{\mu m}``). ```math - PSA = SAN (1-CLA)^{2.4} \\ - PSI = 0.13SIL\\ - PCL = 0.20CLA \\ - SAG = 2.0CLA \, ; \, CLA < 0.25 \\ - SAG = 0.28(CLA-0.25)+0.5 \, ; \, 0.25 \leq CLA \leq 0.5 \\ - SAG = 0.57 \, ; \, CLA > 0.5 \\ - LAG = 1 - PSA - PSI - PCL - SAG + \mathrm{PSA} = \mathrm{SAN} (1-\mathrm{CLA})^{2.4} \\ + \mathrm{PSI} = 0.13\mathrm{SIL}\\ + \mathrm{PCL} = 0.20\mathrm{CLA} \\ + + \mathrm{SAG} = + \begin{align*} + \begin{cases} + 2.0\mathrm{CLA} &\text{ if }\quad \mathrm{CLA} < 0.25 \\ + 0.28(\mathrm{CLA}-0.25)+0.5 &\text{ if }\quad 0.25 \leq \mathrm{CLA} \leq 0.5 \\ + 0.57 &\text{ if }\quad \mathrm{CLA} > 0.5 + \end{cases} + \end{align*} \\ + + \mathrm{LAG} = 1 - \mathrm{PSA} - \mathrm{PSI} - \mathrm{PCL} - \mathrm{SAG} ``` -where ``CLA``, ``SIL`` and ``SAN`` are the primary clay, silt, sand fractions of the topsoil -and ``PCL``, ``PSI``, ``PSA``, ``SAG`` and ``LAG`` are the clay, silt, sand, small and large +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: ```math - TC_{i} = (P_{e})_{i} (S_{g})_{i} \, \rho_{w} \, g \, d_{i} V_{*} + \mathbf{TC}_i = (P_e)_i (S_g)_i \, \rho_w \, g \, d_i V_* ``` -where ``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, ``(S_{g})_{i}`` is the -specific gravity for the particle class i (kg m``^{-3}``), ``\rho_{w}`` is the mass density -of the fluid (kg m``^{-3}``), ``g`` is the acceleration due to gravity (m s``^{-2}``), -``d_{i}`` is the diameter of the particle of class i (m) and ``V_{*}=(g R S)^{0.5}`` is the -shear velocity of the flow (m s``^{-1}``) with ``S`` the slope gradient and ``R`` the -hydraulic radius of the flow (m). The detached sediment are then routed downslope until the -river network using the accucapacityflux, accupacitystate functions depending on the +where ``\mathbf{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 @@ -102,7 +109,7 @@ from land erosion, estimated with the soil loss part of wflow_sediment model, th coming from upstream river cells and the detached sediment that were left in the cell at the end of the previous timestep ``(t-1)``: ```math - (sed_{in})_{t} = (sed_{land})_{t} + upstream\left[(sed_{out})_{t-1}\right] + (sed_{riv})_{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} ``` ### River transport and erosion @@ -126,21 +133,21 @@ Originally more valid for intermediate to large rivers, this simplified version Bagnold equation relates sediment transport to flow velocity with two simple calibration parameters (Neitsch et al, 2011): ```math -C_{max} = c_{sp} \left( \dfrac{prf Q}{h W} \right) ^{sp_{exp}} +C_{\max} = \subtext{c}{sp} \left( \dfrac{\mathrm{prf} Q}{h W} \right)^{\subtext{\mathrm{sp}}{exp}} ``` -where ``C_{max}`` is the sediment concentration (ton m``^{-3}`` or kg/L), ``Q`` is the -surface runoff in the river cell (m``^{3}``s``^{-1}``), ``h`` is the river water level (m), -``W`` is the river width (m) and ``c_{sp}``, ``prf`` and ``sp_{exp}`` are calibration -parameters. The ``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: ``sp_{exp}`` and ``c_{Bagnold} = c_{sp} \, prf^{sp_{exp}}``. The -coefficient ``sp_{exp}`` usually varies between 1 and 2 while ``prf`` and ``c_{sp}`` have a +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) -| Study | River | ``prf`` range | ``c_{sp}`` range | ``sp_{exp}`` range | +| Study | River | ``\mathrm{prf}`` range | ``\subtext{c}{sp}`` range | ``\subtext{\mathrm{sp}}{exp}`` range | |:----- | ----- | ------------- | ---------------- | ------------------ | | Vigiak 2015 | Danube | 0.5-2 (/) | 0.0001-0.01 (0.003-0.006) | 1-2 (1.4) | | Vigiak 2017 | Danube | / | 0.0001-0.01 (0.0015) | 1-2 (1.4) | @@ -151,21 +158,21 @@ Table: Range of the simplified Bagnold coefficients (and calibrated value) models such as Delft3D-WAQ, Engelund and Hansen calculates the total sediment load as (Engelund and Hansen, 1967): ```math - 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} + 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, ``\rho`` and ``\rho_{s}`` are the -fluid and sediment density (here equal to 1000 and 2650 g m``^{-3}``), ``u`` is the water -mean velocity (m s``^{-1}``), ``S`` is the river slope, ``g`` is the acceleration due to gravity, -``D_{50}`` is the river mean diameter (m) and ``\theta`` is the Shields parameter. +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. **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 sediment size. The resulting equation, for a rectangular channel, is (Neitsch et al, 2011): ```math - C_{max} = \left( \dfrac{a u^{b} h^{c} S^{d}}{V_{in}} \right) W + C_{\max} = \left( \dfrac{a u^{b} h^{c} S^{d}}{\subtext{V}{in}} \right) W ``` -where ``V_{in}`` in the volume of water entering the river cell -during the timestep (m``^{3}``) and ``a``, ``b``, ``c`` and ``d`` are coefficients depending +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. @@ -173,28 +180,27 @@ Table: Range of the simplified Bagnold coefficients (and calibrated value) | River sediment diameter | a | b | c | d | |:------------------------|---|---|---|---| -| ``D_{50} \leq`` 0.05mm | 281.4 | 2.622 | 0.182 | 0 | -| 0.05 ``< D_{50} \leq`` 0.25mm | 2 829.6 | 3.646 | 0.406 | 0.412 | -| 0.25 ``< D_{50} \leq`` 2mm | 2 123.4 | 3.300 | 0.468 | 0.613 | -| ``D_{50} >`` 2mm | 431 884.8 | 1.000 | 1.000 | 2.000 | +| ``D_{50} \leq \SI{0.05}{mm}`` | 281.4 | 2.622 | 0.182 | 0 | +| ``\SI{0.05}{mm} < D_{50} \leq \SI{0.25}{mm}`` | 2 829.6 | 3.646 | 0.406 | 0.412 | +| ``\SI{0.25}{mm} < D_{50} \leq \SI{2.0}{mm}`` | 2 123.4 | 3.300 | 0.468 | 0.613 | +| ``D_{50} > \SI{2.0}{mm}`` | 431 884.8 | 1.000 | 1.000 | 2.000 | **Yang** Yang (1996) developed a set of two equations giving transport of sediments for -sand-bed or gravel-bed rivers. The sand equation (``D_{50} < 2mm``) is: +sand-bed or gravel-bed rivers. The sand equation (``D_{50} < \SI{2.0}{mm}``) is: ```math - log\left(C_{ppm}\right) = 5.435 - 0.286log\frac{\omega_{s,50}D_{50}}{\nu}-0.457log\frac{u_{*}}{\omega_{s,50}} \\ - +\left(1.799-0.409log\frac{\omega_{s,50}D_{50}}{\nu}-0.314log\frac{u_{*}}{\omega_{s,50}}\right)log\left(\frac{uS}{\omega_{s,50}}-\frac{u_{cr}S}{\omega_{s,50}}\right) + \log\left(C_{ppm}\right) = 5.435 - 0.286\log\left(\frac{\omega_{s,50}D_{50}}{\nu}\right)-0.457\log\left(\frac{u_*}{\omega_{s,50}}\right) \\ + +\left(1.799-0.409\log\left(\frac{\omega_{s,50}D_{50}}{\nu}\right)-0.314\log\left(\frac{u_*}{\omega_{s,50}}\right)\right)\log\left(\frac{uS}{\omega_{s,50}}-\frac{u_{cr}S}{\omega_{s,50}}\right) ``` -And the gravel equation (``2 \leq D_{50} < 10 mm``) is: +And the gravel equation (``\SI{2.0}{mm} \leq D_{50} < \SI{10.0}{mm}``) is: ```math - log\left(C_{ppm}\right) = 6.681 - 0.633log\frac{\omega_{s,50}D_{50}}{\nu}-4.816log\frac{u_{*}}{\omega_{s,50}} \\ - +\left(2.784-0.305log\frac{\omega_{s,50}D_{50}}{\nu}-0.282log\frac{u_{*}}{\omega_{s,50}}\right)log\left(\frac{uS}{\omega_{s,50}}-\frac{u_{cr}S}{\omega_{s,50}}\right) + \log\left(C_{ppm}\right) = 6.681 - 0.633\log\left(\frac{\omega_{s,50}D_{50}}{\nu}\right)-4.816\log\left(\frac{u_*}{\omega_{s,50}}\right) \\ + +\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, -``\omega_{s,50}`` is the settling velocity of a particle with the median riverbed diameter -estimated with Stokes (m s``^{-1}``), ``\nu`` is the kinematic viscosity of the fluid -(m``^{2}``s``^{-1}``), ``u_{*}`` is the shear velocity (``\sqrt{gR_{H}S}`` in m s``^{-1}`` -with ``R_{H}`` the hydraulic radius of the river) and ``u_{cr}`` is the critical velocity -(m/s, equation can be found in Hessel, 2007). +``\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 sand-bed rivers based on the universal stream power ``\psi``. The corresponding equation is @@ -204,45 +210,45 @@ sand-bed rivers based on the universal stream power ``\psi``. The corresponding ``` where ``\psi`` is the universal stream power given by: ```math - \psi = \dfrac{\psi^{3}}{\left(\dfrac{\rho_{s}}{\rho}-1\right) g h \omega_{s,50} \left[ log_{10}\left(\dfrac{h}{D_{50}}\right)\right]^{2}} + \psi = \dfrac{\psi^{3}}{\left(\dfrac{\rho_{s}}{\rho}-1\right) g h \omega_{s,50} \left[ \log_{10}\left(\dfrac{h}{D_{50}}\right)\right]^{2}} ``` -Once the maximum concentration ``C_{max}`` is established with one of the above transport +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 (``sed_{max} = C_{max} V_{in}``) and the sediment inputs to the river cell -(``sed_{in}`` calculated above) is calculated. If too much sediment is coming in and -``sed_{ex}`` is negative, then there is no river bed and bank erosion. And if the river has +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 ``sed_{stor}`` -are eroded from clay to gravel. If this amount is not enough to cover ``sed_{ex}``, then +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 -``sed_{exeff}`` between ``sed_{ex}`` and ``sed_{stor}``, actual erosion potential is +``\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,bed}`` for the bed and ``E_{R,bank}`` for the river bank. For a +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): ```math - E_{R,bed} = k_{d,bed} \left( \tau_{e,bed} - \tau_{cr,bed} \right) 10^{-6} L W \rho_{b, bed} \Delta t \\~\\ - E_{R,bank} = k_{d,bank} \left( \tau_{e,bank} - \tau_{cr,bank} \right) 10^{-6} L h \rho_{b, bank} \Delta t + E_{R,\mathrm{bed}} = k_{d,\mathrm{bed}} \left( \tau_{e,\mathrm{bed}} - \tau_{cr,\mathrm{bed}} \right) 10^{-6} L W \rho_{b, \mathrm{bed}} \Delta t \\~\\ + 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 ``` -where ``E_{R}`` is the potential bed/bank erosion rates (tons), ``k_{d}`` is the erodibility -of the bed/bank material (cm``^{3}`` N``^{-1}`` s``^{-1}``), ``\tau_{e}`` is the effective -shear stress from the flow on the bed/bank (N m``^{-2}``), ``\tau_{cr}`` is the critical -shear stress for erosion to happen (N m``^{-2}``), ``L``, ``W`` and ``h`` are the channel -length, width and water height (m), ``\rho_{b}`` is the bulk density of the bed/bank of the -river (g cm``^{-3}``) and ``\Delta t`` is the model timestep (s). +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): ```math - k_{d}=0.2 \tau_{cr}^{-0.5} + 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 @@ -251,7 +257,7 @@ efficient enough. The critical shear stress ``\tau_{cr}`` is evaluated different 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 : ```math - \tau_{cr,bank} = (0.1+0.1779 SC+0.0028 SC^{2}-2.34 10^{-5} SC^{3}) C_{ch} + \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 @@ -289,24 +295,23 @@ Then, the repartition of the flow shear stress is refined into the effective she and the bed and bank of the river using the equations developed by Knight (1984) for a rectangular channel: ```math - \tau_{e,bed} = \rho g R_{H} S \left(1 - \dfrac{SF_{bank}}{100}\right) \left(1+\dfrac{2h}{W}\right) \\~\\ - \tau_{e,bank} = \rho g R_{H} S \left( SF_{bank}\right) \left(1+\dfrac{W}{2h}\right) + \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) ``` -where ``\rho g`` is the fluid specific weight (9800 N m``^{-3}`` for water), ``R_{H}`` is the -hydraulic radius of the channel (m), ``h`` and ``W`` are the water level and river width -(m). ``SF_{bank}`` is the proportion of shear stress acting on the bank (%) and is estimated +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): ```math - SF_{bank} = exp \left( -3.230 log_{10}\left(\dfrac{W}{h}+3\right)+6.146 \right) + \mathrm{SF}_{\mathrm{bank}} = \exp \left( -3.230 \log_{10}\left(\dfrac{W}{h}+3\right)+6.146 \right) ``` Finally the relative erosion potential of the bank and bed of the river is calculated by: ```math - RTE_{bed} = \dfrac{E_{R,bed}}{E_{R,bed}+E_{R,bank}} \\~\\ - RTE_{bank} = 1 - RTE_{bed} + \mathrm{RTE}_{\mathrm{bed}} = \dfrac{E_{R,\mathrm{bed}}}{E_{R,\mathrm{bed}}+E_{R,\mathrm{bank}}} \\~\\ + \mathrm{RTE}_{\mathrm{bank}} = 1 - RTE_{\mathrm{bed}} ``` -And the final actual eroded amount for the bed and bank is the maximum between ``RTE -sed_{exeff}`` and the erosion potential ``E_{R}``. Total eroded amount of sediment -``sed_{erod}`` is then the sum of the eroded sediment coming from the storage of previously +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. ### River deposition @@ -315,18 +320,18 @@ the river bed. The deposition process depends on the mass of the sediment, but a characteristics such as velocity. In wflow_sediment, as in SWAT, deposition is modelled with Einstein's equation (Neitsch et al, 2011): ```math - P_{dep}=\left(1-\dfrac{1}{e^{x}}\right)100 + \subtext{P}{dep}=\left(1-\dfrac{1}{e^{x}}\right)100 ``` -where ``P_{dep}`` is the percentage of sediments that is deposited on the river bed and x is +where ``\subtext{P}{dep}`` is the percentage of sediments that is deposited on the river bed and x is a parameter calculated with: ```math x = \dfrac{1.055 L \omega_{s}}{u h} ``` -where ``L`` and ``h`` are channel length and water height (m), ``\omega_{s}`` is the -particle settling velocity calculated with Stokes formula (m s``^{-1}``) and ``u`` is the -mean flow velocity (m s``^{-1}``). The calculated percentage is then subtracted from the -amount of sediment input and eroded river sediment for each particle size class (``sed_{dep} -= P_{dep}/100 (sed_{in} + sed_{erod})``). Resulting deposited sediment are then stored in +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 @@ -334,45 +339,45 @@ Finally after estimating inputs, deposition and erosion with the transport capac flow, the amount of sediment actually leaving the river cell to go downstream is estimated using: ```math - sed_{out} = (sed_{in} + sed_{erod} - sed_{dep}) \dfrac{V_{out}}{V} + \subtext{\mathrm{sed}}{out} = (\subtext{\mathrm{sed}}{in} + \subtext{\mathrm{sed}}{erod} - \subtext{\mathrm{sed}}{dep}) \dfrac{\subtext{V}{out}}{V} ``` -where ``sed_{out}`` is the amount of sediment leaving the river cell (tons), ``sed_{in}`` is +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 in tons), ``sed_{erod}`` is the amount -of sediment coming from river erosion (tons), ``sed_{dep}`` is the amount of deposited -sediments (tons), ``V_{out}`` is the volume of water leaving the river cell (surface runoff -``Q`` times timestep ``\Delta t`` in m``^{3}``) and ``V`` is the total volume of water in -the river cell (``V_{out}`` plus storage ``h W L`` in m``^{3}``). +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 ``(sed_{riv})_{t}``: +end of the timestep ``(\subtext{\mathrm{sed}}{riv})_t``: ```math - (sed_{riv})_{t} = (sed_{riv})_{t-1} + (sed_{land})_{t} + upstream\left[(sed_{out})_{t-1}\right] + (sed_{erod})_{t} - (sed_{dep})_{t} - (sed_{out})_{t} + (\subtext{\mathrm{sed}}{riv})_t = (\subtext{\mathrm{sed}}{riv})_{t-1} + (\subtext{\mathrm{sed}}{land})_t + \mathrm{upstream}\left[(\subtext{\mathrm{sed}}{out})_{t-1}\right] + (\subtext{\mathrm{sed}}{erod})_t - (\subtext{\mathrm{sed}}{dep})_t - (\subtext{\mathrm{sed}}{out})_t ``` ### 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 +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 +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): ```math - TE = \dfrac{\omega_{s}}{u_{cr,res}} = \dfrac{A_{res}}{Q_{out,res}} \omega_{s} + \mathrm{TE} = \dfrac{\omega_s}{u_{cr,\mathrm{res}}} = \dfrac{\subtext{A}{res}}{\subtext{Q}{out,res}} \omega_s ``` -where ``TE`` is the trapping efficiency of the lake/reservoir (or the fraction of particles -trapped), ``\omega_{s}`` is the particle velocity from Stokes (m s``^{-1}``), ``u_{cr,res}`` -is the reservoir's critical settling velocity (m/s) which is equal to the reservoir's -outflow ``Q_{out,res}`` (m``^{3}`` s``^{-1}``) divided by the reservoir's surface area -``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 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: diff --git a/docs/src/model_docs/lateral/waterbodies.md b/docs/src/model_docs/lateral/waterbodies.md index 335ab2867..5b4b4b913 100644 --- a/docs/src/model_docs/lateral/waterbodies.md +++ b/docs/src/model_docs/lateral/waterbodies.md @@ -7,15 +7,14 @@ Simple reservoirs can be included within the river routing by supplying the foll reservoir parameters: + `locs` - Outlet of the reservoirs in which each reservoir has a unique id -+ `area` - Surface area of the reservoirs [m``^2``] ++ `area` - Surface area of the reservoirs ``\SIb{}{m^2}`` + `areas` - Reservoir coverage + `targetfullfrac` - Target fraction full (of max storage) for the reservoir: number between 0 and 1 + `targetminfrac` - Target minimum full fraction (of max storage). Number between 0 and 1 -+ `maxvolume` - Maximum reservoir storage (above which water is spilled) [m``^3``] -+ `demand` - Minimum (environmental) flow requirement downstream of the reservoir [m``^3`` - s``^{-1}``] -+ `maxrelease` - Maximum Q that can be released if below spillway [m``^3`` s``^{-1}``] ++ `maxvolume` - Maximum reservoir storage (above which water is spilled) ``\SIb{}{m^3}`` ++ `demand` - Minimum (environmental) flow requirement downstream of the reservoir ``\SIb{}{m^3 s^{-1}}`` ++ `maxrelease` - Maximum ``Q`` that can be released if below spillway ``\SIb{}{m^3 s^{-1}}`` By default the reservoirs are not included in the model. To include them put the following lines in the TOML file of the model: @@ -44,20 +43,19 @@ targetminfrac = "ResTargetMinFrac" Lakes are modelled using a mass balance approach: ```math - \dfrac{S(t + \Delta t)}{\Delta t} = \dfrac{S(t)}{\Delta t} + Q_{in} + \dfrac{(P-E) A}{\Delta t} - Q_{out} + \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 ``S`` is lake storage [m``^3``], ``\Delta t`` is the model timestep [s], ``Q_{in}`` is -the sum of inflows (river, overland and lateral subsurface flow) [m``^3`` s``^{-1}``], -``Q_{out}`` is the lake outflow at the outlet [m``^3`` s``^{-1}``], ``P`` is precipitation -[m], ``E`` is lake evaporation [m] and ``A`` is the lake surface area [m``^2``]. +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_schematisation](../../images/lake.png) *Lake schematization.* Most of the variables in this equation are already known or coming from previous timestep, -apart from ``S(t+ \Delta t)`` and ``Q_{out}`` which can both be linked to the water level +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 @@ -69,9 +67,9 @@ Approach from Maniak (Burek et al., 2013). Storage curves in wflow can either: 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 ``Q_{out} = \alpha {(H-H_{0})}^{\beta}``, ++ 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 3/2 for a rectangular weir or 2 for a parabolic weir (Bos, 1989). + 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 @@ -79,22 +77,32 @@ relationship between storage and outflow. Storage is assumed to be equal to ``A rating curve for a parabolic weir (``\beta = 2``): ```math - S = A H = A (h + H_{0}) = \dfrac{A}{\sqrt{\alpha}} \sqrt{Q} + A H_{0} + S = A H = A (h + H_{0}) = A \sqrt{\frac{Q}{\alpha}} + A H_0 ``` Inserting this equation in the mass balance gives: ```math - \dfrac{A}{\Delta t \sqrt{\alpha}} \sqrt{Q} + Q = \dfrac{S(t)}{\Delta t} + Q_{in} + - \dfrac{(P-E) A}{\Delta t} - \dfrac{A H_{0}}{\Delta t} = SI - \dfrac{A H_{0}}{\Delta t} + \dfrac{A}{\Delta t} \sqrt{\frac{Q}{\alpha}} + Q = \dfrac{S(t)}{\Delta t} + \subtext{Q}{in} + + A\dfrac{P-E}{\Delta t} - \dfrac{A H_0}{\Delta t} = \mathrm{SI} - \dfrac{A H_0}{\Delta t} ``` -The solution for Q is then: +The solution for ``Q`` is then: ```math - Q = { \left( \dfrac{-LF + \sqrt{LF^{2} + 4 \left( SI - \dfrac{A*H_{0}}{\Delta t} \right)}} - {2} \right) }^{2} \text{for } SI > \dfrac{A H_{0}}{\Delta t} \text{ and where}\\ - LF = \dfrac{A}{\Delta t \sqrt{\alpha}} \\~\\ - Q = 0 \text{ for } SI \leq \dfrac{A*H_{0}}{\Delta t} + 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)} + \right)^2 &\text{ if }\quad \mathrm{SI} > \dfrac{A H_0}{\Delta t} \\ + 0 &\text{ if }\quad \mathrm{SI} \leq \dfrac{A H_0}{\Delta t} + \end{align*} + \end{cases} +``` + +where + +```math + \mathrm{LF} = \dfrac{A}{\Delta t \sqrt{\alpha}}. ``` ### Lake parameters @@ -146,8 +154,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 [m] in the first column `H` and -corresponding lake storage [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 diff --git a/docs/src/model_docs/vertical/sbm.md b/docs/src/model_docs/vertical/sbm.md index 140728566..ca5282c77 100644 --- a/docs/src/model_docs/vertical/sbm.md +++ b/docs/src/model_docs/vertical/sbm.md @@ -1,7 +1,7 @@ # [SBM](@id vert_sbm) ## Introduction -The SBM vertical concept has its roots in the Topog\_SBM model but has had considerable +The SBM vertical concept has its roots in the `Topog_SBM` model but has had considerable changes over time. The main differences are: - The unsaturated zone can be split-up in different layers @@ -42,8 +42,9 @@ 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: + ```math -P'=\frac{-\overline{R}S}{\overline{E}_{w}}ln\left[1-\frac{\overline{E}_{w}}{\overline{R}}(1-p-p_{t})^{-1}\right] +P'=\frac{-\overline{R}S}{\overline{E}_{w}}\log\left[1-\frac{\overline{E}_{w}}{\overline{R}}(1-p-p_{t})^{-1}\right] ``` where ``\overline{R}`` is the average precipitation intensity on a saturated canopy and @@ -105,19 +106,18 @@ path_static = "data/staticmaps-moselle.nc" cyclic = ["vertical.leaf_area_index"] ``` -Furthermore these additional parameters are required: -+ Specific leaf storage (`sl` \[mm\]) -+ Storage woody part of vegetation (`swood` \[mm\]) -+ Extinction coefficient (`kext` \[-\]) +Furthermore, these additional parameters are required: ++ Specific leaf storage (`sl` ``\SIb{}{mm}``) ++ Storage woody part of vegetation (`swood` ``\SIb{}{mm}``) ++ Extinction coefficient (`kext` ``\SIb{}{-}``) -Here it is assumed that `cmax` \[mm\] (leaves) (canopy storage capacity for the leaves only) +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: ```math - - cmax(leaves) = sl \, LAI + c_{\max}(\mathrm{leaves}) = \mathrm{sl} \cdot \mathrm{LAI} ``` To get to total storage (`cmax`) the woody part of the vegetation also needs to be added. As for `sl`, the storage of the woody part `swood` can also be related to land cover (lookup @@ -127,7 +127,7 @@ The canopy gap fraction is determined using the extinction coefficient `kext` (v Bruijnzeel 2001): ```math - canopygapfraction = exp(-kext \, LAI) + \mathrm{canopygapfraction} = \exp(-\subtext{k}{ext} \cdot \mathrm{LAI}) ``` The extinction coefficient `kext` can be related to land cover. @@ -185,18 +185,17 @@ availability. 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` [-]) in each soil layer. A root water uptake reduction model is +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 ``\alpha`` [-] as a function of soil water pressure (``h`` [cm]). Four +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 1 and 5 mm -d``^{-1}``. If ``T_p \le 1 \text{ mm d}^{-1}``, `h3` is set equal to `h3_low` (input model -parameter). If ``T_p \ge 5 \text{ mm d}^{-1}``, `h3` is set equal to `h3_high` (input model -parameter). For ``T_p`` values between 1 and 5 mm d``^{-1}``, the value of `h3` is linearly +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. @@ -205,10 +204,10 @@ The current soil water pressure is determined following the concept defined by B Corey (1964): ```math - \frac{(\theta-\theta_r)}{(\theta_s-\theta_r)} = \Bigg\lbrace{\left(\frac{h_b}{h}\right)^{\lambda}, h > h_b \atop 1 , h \leq h_b} + \frac{\theta-\theta_r}{\theta_s-\theta_r} = \min\left\{1, \left(\frac{h_b}{h}\right)^\lambda\right\} ``` -where ``h`` is the pressure head [cm], ``h_b`` is the air entry pressure head [cm], and +where ``\SIb{h}{cm}`` is the pressure head, ``\SIb{h_b}{cm}`` is the air entry pressure head, and ``\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 @@ -331,35 +330,35 @@ 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 (``z_{t}`` [mm]), divided into a -saturated store (``S`` [mm]) and an unsaturated store (``U`` [mm]). The top of the ``S`` -store forms a pseudo-water table at depth ``z_{i}`` [mm] such that the value of ``S`` at any +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: ```math - S=(z_{t}-z_{i})(\theta_{s}-\theta_{r}) + S=(z_t-z_i)(\theta_s-\theta_r) ``` -where ``\theta_{s}`` [-] and ``\theta_{r}`` [-] are the saturated and residual soil water +where ``\SIb{\theta_{s}}{-}`` and ``\SIb{\theta_{r}}{-}`` are the saturated and residual soil water contents, respectively. -The unsaturated store ``U`` is subdivided into storage (``U_{s}`` [mm]) and deficit -(``U_{d}`` [mm]): +The unsaturated store ``U`` is subdivided into storage (``\SIb{U_s}{mm}``) and deficit +(``\SIb{U_d}{m}``): ```math - U_{d}=(\theta_{s}-\theta_{r})z_{i}-U\\ - U_{s}=U-U_{d} + U_d=(\theta_s-\theta_r)z_i-U\\ + U_s=U-U_d ``` -The saturation deficit (``S_{d}`` [mm]) for the soil profile as a whole is defined as: +The saturation deficit (``\SIb{S_d}{mm}``) for the soil profile as a whole is defined as: ```math - S_{d}=(\theta_{s}-\theta_{r})z_{t}-S + 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 [mm] of the layers in the TOML +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: @@ -368,35 +367,31 @@ mm: thicknesslayers = [100, 300, 800] ``` -The code checks for each grid cell the specified layers against the `soilthickness` [mm], +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 (``st`` [mm t``^{-1}``]) from a ``U`` -[mm] store layer is controlled by the saturated hydraulic conductivity ``K_{sat}`` [mm -t``^{-1}``] at depth ``z`` \[mm\] (bottom layer) or ``z_{i}`` [mm], the effective saturation + +Assuming a unit head gradient, the transfer of water (``\SIb{\mathrm{st}}{mm s^{-1}}``) from a ``\SIb{U}{mm}`` store layer is controlled by the saturated hydraulic conductivity ``\SIb{\subtext{K}{sat}}{mm s^{-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): ```math - st=K_{\mathit{sat}}\left(\frac{\theta-\theta_{r}}{\theta_{s}-\theta_{r}}\right)^{c}\\~\\ + \mathrm{st}=\subtext{K}{sat}\left(\frac{\theta-\theta_r}{\theta_s-\theta_r}\right)^c\\~\\ c=\frac{2+3\lambda}{\lambda} ``` When the unsaturated layer is not split-up into different layers, it is possible to use the -original Topog\_SBM vertical transfer formulation, by specifying in the TOML file: +original `Topog\_SBM` vertical transfer formulation, by specifying in the TOML file: ```toml [model] transfermethod = true ``` -The transfer of water from the ``U`` [mm] store to the ``S`` [mm] store (``st`` [mm -t``^{-1}``]) is in that case controlled by the saturated hydraulic conductivity ``K_{sat}`` -[mm t``^{-1}``] at depth ``z_{i}`` [mm] and the ratio between ``U`` [mm] and ``S_{d}`` -[mm]: +The transfer of water from the ``\SIb{U}{mm}`` store to the ``\SIb{S}{mm}`` store (``\SIb{st}{mm s^{-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}``: ```math - st=K_{\mathit{sat}}\frac{U_{s}}{S_{d}} + \mathrm{st}=\subtext{K}{sat}\frac{U_s}{S_d} ``` Four different saturated hydraulic conductivity depth profiles (`ksat_profile`) are @@ -407,65 +402,65 @@ available and a `ksat_profile` can be specified in the TOML file as follows: ksat_profile = "exponential_constant" # optional, one of ("exponential", "exponential_constant", "layered", "layered_exponential"), default is "exponential" ``` -Soil measurements are often available for about the upper 1.5-2 m of the soil column to +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 1.5-2 m. These different profiles allow to extent the saturated +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 (``K_{sat}`` [mm -t``^{-1}``]) declines with soil depth (``z`` [mm]) in the model according to: +`ksat_profile` "exponential", the saturated hydraulic conductivity ``\SIb{\subtext{K}{sat}}{mm s^{-1}}`` declines with soil depth ``\SIb{z}{mm}`` in the model according to: ```math - K_{sat}=K_{0}e^{(-fz)}, + \subtext{K}{sat} = K_0 e^{-fz}, ``` -where ``K_{0}`` [mm t``^{-1}``] is the saturated hydraulic conductivity at the soil surface -and ``f`` is a scaling parameter [mm``^{-1}``]. +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 ``K_{sat}`` for different values of ``f``. +conductivity ``\subtext{K}{sat}`` for different values of ``f``. ```@setup plot using Printf using CairoMakie ``` + ```@example plot let # hide fig = Figure(resolution = (800, 400)) # hide - ax = Axis(fig[1, 1], xlabel = "Kₛₐₜ [mm/day]", ylabel = "-z [mm]") # hide + ax = Axis(fig[1, 1], xlabel = L"K_\mathrm{sat}\;[\mathrm{mm/day}]", ylabel = L"-z\;[\mathrm{mm}]") # hide z = 0:5.0:1000 # hide ksat = 100.0 # hide f = 0.6 ./ collect(50:150.0:800) # hide for fi in f # hide - lines!(ax, ksat .* exp.(-fi .* z), -z, label = @sprintf("f = %.2e", fi)) # hide + lines!(ax, ksat .* exp.(-fi .* z), -z, label = @sprintf("%.2e", fi)) # hide end # hide - Legend(fig[1, 2], ax, "f") # hide + Legend(fig[1, 2], ax, L"f") # hide fig # hide end # hide ``` -With `ksat_profile` "exponential\_constant", ``K_{sat}`` declines exponentially with soil -depth ``z`` until ``z_\mathrm{exp}`` [mm] below the soil surface, and stays constant at and -beyond soil depth ``z_\mathrm{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}``: ```math - K_{sat} = \begin{cases} - K_{0}e^{(-fz)} & \text{if $z < z_\mathrm{exp}$}\\ - K_{0}e^{(-fz_\mathrm{exp})} & \text{if $z \ge z_\mathrm{exp}$}. + \subtext{K}{sat} = \begin{cases} + K_0e^{-fz} & \text{if $z < \subtext{z}{exp}$}\\ + K_0e^{-f\subtext{z}{exp}} & \text{if $z \ge \subtext{z}{exp}$}. \end{cases} ``` -It is also possible to provide a ``K_{sat}`` value per soil layer by specifying -`ksat_profile` "layered", these ``K_{sat}`` values are used directly to compute the vertical +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 ``K_{sat}`` value per soil layer is used until depth -``z_\mathrm{layered}`` below the soil surface, and beyond ``z_\mathrm{layered}`` an -exponential decline of ``K_{sat}`` (of the soil layer with bottom ``z_\mathrm{layered}``) +`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 ``K_{sat}`` value. +the figure below where the blue line represents the ``\subtext{K}{sat}`` value. ![ksat_profiles](../../images/sbm_ksat_profiles.png) @@ -479,6 +474,7 @@ have different infiltration capacities. Naturally, only the water that can be st 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 @@ -507,13 +503,12 @@ The near surface soil temperature is modelled using a simple equation (Wigmosta 2009): ```math -T_s^{t} = T_s^{t-1} + w (T_a - T_s^{t-1}) +T_s^t = T_s^{t-1} + w (T_a - T_s^{t-1}) ``` -where ``T_s^{t}`` [``\degree``C] is the near-surface soil temperature at time ``t``, ``T_a`` -[``\degree``C] is air temperature and ``w`` [-] is a weighting coefficient determined +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` [-], default is 0.038) is applied to the maximum infiltration +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: @@ -530,8 +525,8 @@ A S-curve (see plot below) is used to make a smooth transition (a c-factor (``c` used): ```math - b = \frac{1.0}{(1.0 - cf\_soil)}\\~\\ - soilinfredu = \frac{1.0}{b + exp(-c (T_s - a))} + cf\_soil\\~\\ + b = \frac{1.0}{1.0 - \subtext{\mathrm{cf}}{soil}}\\~\\ + \mathrm{soilinfredu} = \frac{1.0}{b + \exp(-c (T_s - a))} + \subtext{\mathrm{cf}}{soil}\\~\\ a = 0.0\\ c = 8.0 ``` @@ -548,12 +543,13 @@ used): ### Capillary rise + The actual capillary rise `actcapflux` [mm t``^{-1}``] is determined using the following approach: first the saturated hydraulic conductivity `ksat` [mm t``^{-1}``] is determined at -the water table ``z_{i}``; next a potential capillary rise `maxcapflux` [mm t``^{-1}``] is +the water table ``z_i``; next a potential capillary rise `maxcapflux` [mm t``^{-1}``] is determined from the minimum of `ksat`, actual transpiration `actevapustore` [mm t``^{-1}``] -taken from the ``U`` store, available water in the ``S`` store (`satwaterdepth` [mm]) and -the deficit of the ``U`` store (`ustorecapacity` [mm]), as shown by the following code +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 block: ```julia @@ -561,8 +557,7 @@ block: ``` Then the potential rise `maxcapflux` is scaled using the water table depth `zi`, a maximum -water depth `cap_hmax` [mm] beyond which capillary rise ceases and a coefficient `cap_n` -[-], as follows in the code block below (`i` refers to the index of the vector that contains +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 @@ -589,8 +584,8 @@ the layer position): usl[k] * (sbm.theta_s[i] - sbm.theta_r[i]) - usld[k] ``` -where `usl` [mm] is the unsaturated layer thickness, `usld` is the `ustorelayerdepth` \[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 @@ -615,7 +610,7 @@ capillary rise `netcapflux` [mm t``^{-1}``]. ### Leakage -If the `maxleakage` (mm/day) input model parameter is set > 0, water is lost from the +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 @@ -645,7 +640,7 @@ returned as surface water. The return flow fraction (``f_\mathrm{return}`` [-]) calculated as follows: ```math - f_\mathrm{return} = 1.0 - \frac{d_\mathrm{net}}{d_\mathrm{gross}}, + \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 @@ -661,21 +656,21 @@ supported. These computations can be enabled by specifying the following in the nonpaddy = true ``` Irrigation is applied during the growing season (when input parameter `irrigation_trigger` -[-] is `true` (or `on`)) and when water depletion exceeds the readily available water: +``\SIb{}{-}`` is `true` (or `on`)) and when water depletion exceeds the readily available water: ```math - (U_\mathrm{field} - U_\mathrm{a}) \ge (U_\mathrm{field} - U_\mathrm{h3}) + (\subtext{U}{field} - \subtext{U}{a}) \ge (\subtext{U}{field} - \subtext{U}{h3}) ``` -where ``U_\mathrm{field}`` \[mm\] is the unsaturated store in the root zone at field -capacity (defined at a soil water pressure head of -100 cm), ``U_\mathrm{a}`` \[mm\] is the -actual unsaturated store in the root zone and ``U_\mathrm{h3}`` \[mm\] is the unsaturated +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 [mm t``^{-1}``] is the irrigation rate that brings the root zone back to field capacity, limited by the soil infiltration capacity [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` [-], default is 1.0), resulting in gross irrigation +non-paddy crops (`irrigation_efficiency` ``\SIb{}{-}}``, default is ``1.0``), resulting in gross irrigation demand [mm t``^{-1}``]. Finally, the gross irrigation demand is limited by the maximum irrigation rate (`maximum_irrigation_rate` [mm t``^{-1}``], default is 25 mm d``^{-1}``). If the maximum irrigation rate is applied, irrigation continues at subsequent time steps until @@ -691,26 +686,25 @@ computations can be enabled by specifying the following in the TOML file: paddy = true ``` Irrigation is applied during the growing season (when input parameter `irrigation_trigger` -[-] is `true` (or `on`)) and when the paddy water depth `h` \[mm\] reaches below the minimum -water depth `h_min` \[mm\] (see also the figure below). The net irrigation demand [mm +``\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` -\[mm\], an approach similar to Xie and Cui (2011). To account for limited irrigation +``\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` [-], default is 1.0), resulting in gross irrigation demand +fields (`irrigation_efficiency` ``\SIb{}{-}``, default is 1.0), resulting in gross irrigation demand [mm t``^{-1}``]. Finally, the gross irrigation demand is limited by the maximum irrigation -rate (`maximum_irrigation_rate` [mm t``^{-1}``], default is 25 mm d``^{-1}``). If the +rate (`maximum_irrigation_rate` [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` [mm t``^{-1}``], the amount of water available for infiltration. When the -paddy water depth `h` exceeds `h_max` \[mm\] runoff occurs, and this amount is added to 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 5 mm d``^{-1}``), and a +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` [-], a -multiplication factor applied to the vertical hydraulic conductivity at soil depth ``z`` -[mm]. +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) @@ -734,61 +728,57 @@ 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 +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` -[-], described in the next sub-section. +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 ``V_\mathrm{sw, demand}`` [m``^3``] and availability -``V_\mathrm{sw, availabilty}`` [m``^3``] are summed (including reservoir and lake locations -limited by a fixed scaling factor of 0.98), and the total surface water abstraction is then: +For allocation areas the water demand ``\SIb{\subtext{V}{sw, demand}}{m^3}`` and availability +``\SIb{\subtext{V}{sw, availabilty}}{m^3}`` are summed (including reservoir and lake locations +limited by a fixed scaling factor of ``0.98``), and the total surface water abstraction is then: ```math - V_\mathrm{sw, abstraction} = \mathrm{min}(V_\mathrm{sw, demand}, V_\mathrm{sw, availabilty}) + \subtext{V}{sw, abstraction} = \min (\subtext{V}{sw, demand}, \subtext{V}{sw, availabilty}) ``` -The fraction of available surface water that can be abstracted ``f_\mathrm{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: ```math - f_\mathrm{sw, abstraction} = \frac{V_\mathrm{sw, abstraction}}{V_\mathrm{sw, available}} + \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. The fraction of water demand that can be satisfied by available surface water -``f_\mathrm{sw, allocation}`` [-] at the allocation area level is then: +``\SIb{\subtext{f}{sw, allocation}}{-}`` at the allocation area level is then: ```math - f_\mathrm{sw, allocation} = \frac{V_\mathrm{sw, abstraction}}{V_\mathrm{sw, demand}} + \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. 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 ``V_\mathrm{gw, -demand}`` [m``^3``] and availability ``V_\mathrm{gw, availabilty}`` [m``^3``] are summed, +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: ```math - V_\mathrm{gw, abstraction} = \mathrm{min}(V_\mathrm{gw, demand}, V_\mathrm{gw, availabilty}) + \subtext{V}{gw, abstraction} = \min(\subtext{V}{gw, demand}, \subtext{V}{gw, availabilty}) ``` The fraction of available groundwater that can be abstracted at allocation area level -``f_\mathrm{gw, abstraction}`` [-] at the allocation area level is then: +``\SIb{\subtext{f}{gw, abstraction}}{-}`` at the allocation area level is then: ```math - f_\mathrm{gw, abstraction} = \frac{V_\mathrm{gw, abstraction}}{V_\mathrm{gw, available}} + \subtext{f}{gw, abstraction} = \frac{\subtext{V}{gw, abstraction}}{\subtext{V}{gw, available}} ``` 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 ``f_\mathrm{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: ```math - f_\mathrm{gw, allocation} = \frac{V_\mathrm{gw, abstraction}}{V_\mathrm{gw, demand}} + \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. @@ -796,14 +786,14 @@ 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` [m``^3`` s``^{-1}``] is divided by -the flow length `dl` [m] and subtracted from the lateral inflow of kinematic wave routing +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` [m``^3`` s``^{-1}``] is subtracted as part of the +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` [m``^3`` t``^{-1}``]) from the reservoir -`volume` [m``^3``] and lake `storage` [m``^3``] respectively, with a subsequent update of -the lake `waterlevel` [m]. +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