Skip to content

Commit

Permalink
Support different vertical saturated hydraulic conductivity profiles …
Browse files Browse the repository at this point in the history 
…in `SBM` (#340)

* Ksat soil profiles initialization
Exponential and exponential-constant profiles, based on `z_exp` (depth for which exponential decline of ksat is valid) and profile setting in TOML file (default is `exponential`).

* Ksat profile vertical flow

Co-authored-by: Raul Mendoza <[email protected]>

* Ksat profiles for lateral subsurface flow
Exponential and exponential-constant profile.

Co-authored-by: Raul Mendoza <[email protected]>

* Additonal ksat profiles vertical flow
`layered` and `layered_exponential`

Co-authored-by: Raul Mendoza <[email protected]>

* Additonal ksat profiles lateral subsurface flow
Initialization of `layered` and `layered_exponential` profiles for lateral subsurface flow and computation of equivalent horizontal conductivity for layered profile.

Co-authored-by: Raul Mendoza <[email protected]>

* Kinematic wave ssf and ksat profiles
Add kinematic wave ssf for ksat profile `layered` and `layered_exponential`.

Co-authored-by: Raul Mendoza <[email protected]>

* Update capflux and deeptransfer for ksatprofiles
Computation of ksat at zi (capflux) and deepksat at bottom soil (deeptransfer).

* Updates and fixes related to `ksat_profile`

* Add tests and couple of fixes
Fixes are related to ksat_profile (vertical and horizontal hydraulic conductivity).

* Number of layers ksat_profile
For `layered` profile use `nlayers` for number of layers with `kv` values, and for `layered` and `layered_exponential` profile use `nlayers` to compute `ssfmax`.

* Update docs `SBM`

* Update docs lateral subsurface flow

* Update changelog
And fix small typo docs.

* Refactor initialization of lateral subsurface flow based on `ksat_profile`
Make use of smaller subfunctions.

* Add ksat_profile options to docs lateral subsurface flow

* Add variable `z_layered` for `ksat_profile`
For the layered_exponential profile `z_exp` was used, this is the depth from the soil surface for which an exponential decline is valid. For the layered_exponential profile the use of `z_exp` is confusing and  `z_layered` (depth from soil surface for which a layered profile is valid) is a better term.

* Add `z_layered` to params_vertical.md

* Add `ksat_profile` conceptual figure to docs

---------

Co-authored-by: Raul Mendoza <[email protected]>
  • Loading branch information
@vers-w @raulmendoza10
vers-w and raulmendoza10 authored Mar 12, 2024
1 parent 27c3160 commit 7be5558
Show file tree
Hide file tree
Showing 14 changed files with 699 additions and 76 deletions.
3 changes: 3 additions & 0 deletions docs/src/changelog.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- Checks to see if all states are covered in the .toml file. If not all states are covered,
an error is thrown. If there are more states specified than required, these states are
ignored (with a warning in the logging), and the simulation will continue.
- Support different vertical hydraulic conductivity profiles for the `SBM` concept. See also
the following sections: [The SBM soil water accounting scheme](@ref) and [Subsurface flow
routing](@ref) for a short description.

## v0.7.3 - 2024-01-12

Expand Down
Binary file added docs/src/images/sbm_ksat_profiles.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
64 changes: 49 additions & 15 deletions docs/src/model_docs/lateral/kinwave.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,30 +55,64 @@ kinematic wave for surface water routing, as a cyclic parameter or as part of fo
also [Input section](@ref)).

## Subsurface flow routing
In the SBM model the kinematic wave approach is used to route subsurface flow laterally. The
saturated store ``S`` can be drained laterally by saturated downslope subsurface flow per
unit width of slope ``w`` [m] according to:
In the SBM model the kinematic wave approach is used to route subsurface flow laterally.
Different vertical hydraulic conductivity depth profiles are possible as part of the
vertical [SBM](@ref soil) concept, and these profiles (after unit conversion) are also used
to compute lateral subsurface flow. The following profiles (see [SBM](@ref soil) for a
detailed description) are available:
- `exponential` (default)
- `exponential_constant`
- `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]
according to:
```math
q=\frac{K_{0}\mathit{tan(\beta)}}{f}(e^{(-fz_{i})}-e^{(-fz_{t})})
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}$}\\
\\
K_0e^{(-fz_\mathrm{exp})}(z_t - z_i)\tan(\beta) w & \text{if $z_i \ge z_\mathrm{exp}$},
\end{cases}
```
where ``\beta`` is element slope angle [deg.], ``q`` is subsurface flow [m``^{2}``/t],
``K_{0}`` is the saturated hydraulic conductivity at the soil surface [m/t], ``z_{i}`` is
the water table depth [m], ``z_{t}`` is total soil depth [m], and ``f`` is a scaling
parameter [m``^{-1}``], that controls the decrease of vertical saturated conductivity with
depth.
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``.

Combining with the following continuity equation:
```math
(\theta_s-\theta_r)\frac{\partial h}{\partial t} = -w\frac{\partial q}{\partial x} + wr
(\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/t] to the saturated store. Substituting for ``h (\frac{\partial
q}{\partial h})``, gives:
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
(\frac{\partial Q}{\partial h})``, gives:
```math
w \frac{\partial q}{\partial t} = -cw\frac{\partial q}{\partial x} + cwr
\frac{\partial Q}{\partial t} = -c\frac{\partial Q}{\partial x} + cwr
```

where celerity ``c = \frac{K_{0}\mathit{tan(\beta)}}{(\theta_s-\theta_r)} e^{(-fz_{i})}``
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}$}\\
\\
\frac{K_0e^{(-fz_\mathrm{exp})}\tan(\beta)}{(\theta_s-\theta_r)} & \text{if $z_i \ge z_\mathrm{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:
```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)}.
```

The kinematic wave equation for lateral subsurface flow is solved iteratively using Newton's
method.
Expand Down
31 changes: 21 additions & 10 deletions docs/src/model_docs/params_lateral.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,11 +155,11 @@ between parentheses.
The Table below shows the parameters (fields) of struct `LateralSSF`, including a
description of these parameters, the unit, and default value if applicable. The parameters
in bold represent model parameters that can be set through static input data (netCDF). The
soil related parameters `f`, `soilthickness`, `θₛ` and `θᵣ` are derived from the vertical
`SBM` concept (including unit conversion for `f` and `soilthickness`), and can be listed in
the TOML configuration file under `[input.vertical]`, to map the internal model parameter to
the external netCDF variable. The internal slope model parameter `βₗ` is set through the
TOML file as follows:
soil related parameters `f`, `soilthickness`, `z_exp`, `θₛ` and `θᵣ` are derived from the
vertical `SBM` concept (including unit conversion for `f`, `z_exp` and `soilthickness`), and
can be listed in the TOML configuration file under `[input.vertical]`, to map the internal
model parameter to the external netCDF variable. The internal slope model parameter `βₗ` is
set through the TOML file as follows:

```toml
[input.lateral.land]
Expand All @@ -168,22 +168,33 @@ slope = "Slope"

The parameter `kh₀` is computed by multiplying the vertical hydraulic conductivity at the
soil surface `kv₀` (including unit conversion) of the vertical `SBM` concept with the
external parameter `ksathorfrac` \[-\] (default value of 1.0). The external parameter
`ksathorfrac` can be set as follows through the TOML file:
internal parameter `khfrac` \[-\] (default value of 1.0). The internal model parameter
`khfrac` is set through the TOML file as follows:

```toml
[input.lateral.subsurface]
ksathorfrac = "KsatHorFrac"
```

The `ksathorfrac` parameter compensates for anisotropy, small scale `kv₀` measurements (soil
The `khfrac` parameter compensates for anisotropy, small scale `kv₀` measurements (soil
core) that do not represent larger scale hydraulic conductivity, and smaller flow length
scales (hillslope) in reality, not represented by the model resolution.

For the vertical [SBM](@ref params_sbm) concept different vertical hydraulic conductivity
depth profiles are possible, and these also determine which `LateralSSF` parameters are used
including the input requirements for the computation of lateral subsurface flow. For the
`exponential` profile the model parameters `kh₀` and `f` are used. For the
`exponential_constant` profile `kh₀` and `f` are used, and `z_exp` is required as part of
`[input.vertical]`. For the `layered` profile, `SBM` model parameter `kv` is used, and for
the `layered_exponential` profile `kv` is used and `z_exp` is required as part of
`[input.vertical]`.

| parameter | description | unit | default |
|:---------------| --------------- | ---------------------- | ----- |
| `kh₀` | horizontal hydraulic conductivity at soil surface | m d``^{-1}`` | 3.0 |
| **`f`** | a scaling parameter (controls exponential decline of kh₀) | m``^{-1}`` | 1.0 |
| **`f`** | a scaling parameter (controls exponential decline of `kh₀`) | m``^{-1}`` | 1.0 |
| `kh` | horizontal hydraulic conductivity | m d``^{-1}`` | - |
| **`khfrac`** (`ksathorfrac`) | a muliplication factor applied to vertical hydraulic conductivity `kv` | - | 100.0 |
| **`soilthickness`** | soil thickness | m | 2.0 |
| **`θₛ`** | saturated water content (porosity) | - | 0.6 |
| **`θᵣ`** | residual water content | - | 0.01 |
Expand All @@ -192,13 +203,13 @@ scales (hillslope) in reality, not represented by the model resolution.
| `dl` | drain length | m | - |
| `dw` | drain width | m | - |
| `zi` | pseudo-water table depth (top of the saturated zone) | m | - |
| **`z_exp`** | depth from soil surface for which exponential decline of `kh₀` is valid | m | - |
| `exfiltwater` | exfiltration (groundwater above surface level, saturated excess conditions) | m Δt⁻¹ | - |
| `recharge` | net recharge to saturated store | m``^2`` Δt⁻¹ | - |
| `ssf` | subsurface flow | m``^3`` d``{-1}`` | - |
| `ssfin` | inflow from upstream cells | m``^3`` d``{-1}`` | - |
| `ssfmax` | maximum subsurface flow | m``^2`` d``{-1}`` | - |
| `to_river` | part of subsurface flow that flows to the river | m``^3`` d``{-1}`` | - |
| `wb_pit` | boolean location (0 or 1) of a waterbody (wb, reservoir or lake) | - | - |

## Local inertial

Expand Down
22 changes: 19 additions & 3 deletions docs/src/model_docs/params_vertical.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,20 @@ external netCDF variable `Sl`:
specific_leaf = "Sl"
```

Different [vertical hydraulic conductivity depth profiles](@ref soil): `exponential`
(default), `exponential_constant`, `layered` and `layered_exponential` can be provided
through the TOML file. Below an example for the `exponential_constant` profile:

```toml
[input.vertical]
ksat_profile = "exponential_constant"
```

For the `exponential` profile the input parameters `kv₀` and `f` are used. For the
`exponential_constant` profile `kv₀` and `f` are used, and `z_exp` is required as input. For
the `layered` profile, input parameter `kv` is used, and for the `layered_exponential`
profile `kv` is used and `z_layered` is required as input.

| parameter | description | unit | default |
|:---------------| --------------- | ---------------------- | ----- |
| **`cfmax`** | degree-day factor | mm ᵒC``^{-1}`` Δt``^{-1}`` | 3.75653 mm ᵒC``^{-1}`` day``^{-1}`` |
Expand All @@ -33,7 +47,10 @@ specific_leaf = "Sl"
| **`θₛ`** (`theta_s`) | saturated water content (porosity) | - | 0.6 |
| **`θᵣ`** (`theta_r`) | residual water content | - | 0.01 |
| **`kv₀`** (`kv_0`) | Vertical hydraulic conductivity at soil surface | mm Δt``^{-1}`` | 3000.0 mm day``^{-1}``|
| **`f`** | scaling parameter (controls exponential decline of kv₀) | mm``^{-1}`` | 0.001 |
| **`kv`** | Vertical hydraulic conductivity per soil layer | mm Δt``^{-1}`` | 1000.0 mm day``^{-1}``|
| **`f`** | scaling parameter (controls exponential decline of `kv₀`) | mm``^{-1}`` | 0.001 |
| **`z_exp`** | Depth from soil surface for which exponential decline of `kv₀` is valid | mm | - |
| **`z_layered`** | Depth from soil surface for which layered profile (of `layered_exponential`) is valid | mm | - |
| **`hb`** | air entry pressure of soil (Brooks-Corey) | cm | 10.0 |
| **`soilthickness`** | soil thickness | mm | 2000.0 |
| **`infiltcappath`** | infiltration capacity of the compacted areas | mm Δt``^{-1}`` | 10.0 mm day``^{-1}`` |
Expand All @@ -59,6 +76,7 @@ specific_leaf = "Sl"
| `n` | number of grid cells | - | - |
| `nlayers` | number of soil layers | - | - |
| `n_unsatlayers` | number of unsaturated soil layers | - | - |
| `nlayers_kv` | number of soil layers with vertical hydraulic conductivity value `kv` | - | - |
| `riverfrac` | fraction of river | - | - |
| `act_thickl` | thickness of soil layers | mm | - |
| `sumlayers` | cumulative sum of soil layers thickness, starting at soil surface | mm | - |
Expand All @@ -69,7 +87,6 @@ specific_leaf = "Sl"
| `zi`| pseudo-water table depth (top of the saturated zone) | mm | - |
| `soilwatercapacity`| soilwater capacity | mm | - |
| `canopystorage`| canopy storage | mm | - |
| `canopygapfraction` | canopygapfraction | - | - |
|**`precipitation`** | precipitation | mm Δt``^{-1}``| - |
| **`temperature`** | temperature | ᵒC | - |
| **`potential_evaporation`** | potential evaporation | mm Δt``^{-1}`` | - |
Expand Down Expand Up @@ -111,7 +128,6 @@ specific_leaf = "Sl"
| `snow` | snow storage | mm | - |
| `snowwater` | liquid water content in the snow pack | mm | - |
| `rainfallplusmelt` | snowmelt + precipitation as rainfall | mm Δt``^{-1}`` | - |
| `glacierstore` | water within the glacier | mm | - |
| `tsoil` | top soil temperature | ᵒC | - |
| **`leaf_area_index`** | leaf area index | m``^2`` m``{-2}`` | - |
| `waterlevel_land` | water level land | mm | - |
Expand Down
48 changes: 43 additions & 5 deletions docs/src/model_docs/vertical/sbm.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -350,14 +350,28 @@ t``^{-1}``]) is in that case controlled by the saturated hydraulic conductivity
```math
st=K_{\mathit{sat}}\frac{U_{s}}{S_{d}}
```
Saturated conductivity (``K_{sat}`` [mm t``^{-1}``]) declines with soil depth (``z`` [mm])
in the model according to:

Four different saturated hydraulic conductivity depth profiles (`ksat_profile`) are
available and a `ksat_profile` can be specified in the TOML file as follows:

```toml
[input.vertical]
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
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
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:

```math
K_{sat}=K_{0}e^{(-fz)}
K_{sat}=K_{0}e^{(-fz)},
```
where ``K_{0}`` [mm t``^{-1}``] is the saturated conductivity at the soil surface and ``f``
is a scaling parameter [mm``^{-1}``].
where ``K_{0}`` [mm t``^{-1}``] is the saturated hydraulic conductivity at the soil surface
and ``f`` is a scaling parameter [mm``^{-1}``].

The plot below shows the relation between soil depth ``z`` and saturated hydraulic
conductivity ``K_{sat}`` for different values of ``f``.
Expand Down Expand Up @@ -385,6 +399,30 @@ conductivity ``K_{sat}`` for different values of ``f``.
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}``:

```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}$}.
\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
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}``)
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.

![ksat_profiles](../../images/sbm_ksat_profiles.png)

*Overview of available `ksat_profile` options, for a soil column with five layers*

### Infiltration

The water available for infiltration is taken as the rainfall including meltwater.
Expand Down
57 changes: 40 additions & 17 deletions src/flow.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -389,6 +389,8 @@ end
@get_units @exchange @grid_type @grid_location @with_kw struct LateralSSF{T}
kh₀::Vector{T} | "m d-1" # Horizontal hydraulic conductivity at soil surface [m d⁻¹]
f::Vector{T} | "m-1" # A scaling parameter [m⁻¹] (controls exponential decline of kh₀)
kh::Vector{T} | "m d-1" # Horizontal hydraulic conductivity [m d⁻¹]
khfrac::Vector{T} | "-" # A muliplication factor applied to vertical hydraulic conductivity `kv` [-]
soilthickness::Vector{T} | "m" # Soil thickness [m]
θₛ::Vector{T} | "-" # Saturated water content (porosity) [-]
θᵣ::Vector{T} | "-" # Residual water content [-]
Expand All @@ -397,6 +399,7 @@ end
dl::Vector{T} | "m" # Drain length [m]
dw::Vector{T} | "m" # Flow width [m]
zi::Vector{T} | "m" # Pseudo-water table depth [m] (top of the saturated zone)
z_exp::Vector{T} | "m" # Depth [m] from soil surface for which exponential decline of kv₀ is valid
exfiltwater::Vector{T} | "m Δt-1" # Exfiltration [m Δt⁻¹] (groundwater above surface level, saturated excess conditions)
recharge::Vector{T} | "m2 Δt-1" # Net recharge to saturated store [m² Δt⁻¹]
ssf::Vector{T} | "m3 d-1" # Subsurface flow [m³ d⁻¹]
Expand All @@ -411,7 +414,7 @@ end
end


function update(ssf::LateralSSF, network, frac_toriver)
function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
@unpack subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes = network

ns = length(subdomain_order)
Expand All @@ -433,22 +436,42 @@ function update(ssf::LateralSSF, network, frac_toriver)
upstream_nodes[n],
eltype(ssf.to_river),
)

ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf(
ssf.ssfin[v],
ssf.ssf[v],
ssf.zi[v],
ssf.recharge[v],
ssf.kh₀[v],
ssf.βₗ[v],
ssf.θₛ[v] - ssf.θᵣ[v],
ssf.f[v],
ssf.soilthickness[v],
ssf.Δt,
ssf.dl[v],
ssf.dw[v],
ssf.ssfmax[v],
)
if (ksat_profile == "exponential") ||
(ksat_profile == "exponential_constant")
ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf(
ssf.ssfin[v],
ssf.ssf[v],
ssf.zi[v],
ssf.recharge[v],
ssf.kh₀[v],
ssf.βₗ[v],
ssf.θₛ[v] - ssf.θᵣ[v],
ssf.f[v],
ssf.soilthickness[v],
ssf.Δt,
ssf.dl[v],
ssf.dw[v],
ssf.ssfmax[v],
ssf.z_exp[v],
ksat_profile,
)
elseif (ksat_profile == "layered") ||
(ksat_profile == "layered_exponential")
ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf(
ssf.ssfin[v],
ssf.ssf[v],
ssf.zi[v],
ssf.recharge[v],
ssf.kh[v],
ssf.βₗ[v],
ssf.θₛ[v] - ssf.θᵣ[v],
ssf.soilthickness[v],
ssf.Δt,
ssf.dl[v],
ssf.dw[v],
ssf.ssfmax[v],
)
end
end
end
end
Expand Down
Loading

0 comments on commit 7be5558

Please sign in to comment.