Merge branch 'master' into refactor/sediment_concept

docs/changelog.qmd

@@ -16,6 +16,11 @@ project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 - Mutating BMI model control functions (`update`, `update_until` and `finalize`) and
   extended mutating BMI functions (`load_state` and `save_state`) should return `nothing`.
 - Added downloading of testdata to Dockerfile, to ensure an image was able to build.
+- The reservoir (`reservoir_index_f`) and lake (`lake_index_f`) indices as part of
+  `network.river` were not correct. These were mapped to their own index in the
+  `SimpleReservoir` and `Lake` struct, and not to the corresponding river index. This
+  resulted in incorrect surface water abstractions from reservoir and lake volumes, and
+  surface water abstractions were set at zero at the wrong river locations.
 
 ### Changed
 - Removed vertical concepts `HBV` and `FLEXTopo`.
@@ -31,6 +36,13 @@ project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   model `variables`, `parameters` and `boundary_conditions` (if applicable), including
   associated functions for initializing and updating these model components. The original
   long update function of the `SBM` soil part has been split into separate functions.
+- Refactor the lateral (routing) components: as for the vertical `SBM` concept split the
+  structs into `variables`, `parameters` and `boundary_conditions` (if applicable).
+- Timestepping method parameters for solving the kinematic wave and local inertial
+  approaches for river and overland flow are moved to a `TimeStepping` struct. The
+  timestepping implementation for the kinematic wave is now similar to the local inertial
+  method: a stable timestep is computed for each sub timestep (or a fixed sub timestep is
+  used) as part of a while loop (for each model timestep).
 
 ### Added
 - Support direct output of snow and glacier melt, and add computation of snow water

docs/home/publications.qmd

@@ -148,6 +148,11 @@ tracer data: Application to the data-scarce area of the Bandung groundwater basi
 Indonesia, Journal of Hydrology: Regional Studies, 50,
 <https://doi.org/10.1016/j.ejrh.2023.101585>.
 
+Rusli, S. R., Bense, V. F., Mustafa, S. M. T., and Weerts, A. H.,2024. The impact of future
+changes in climate variables and groundwater abstraction on basin-scale groundwater
+availability, Hydrol. Earth Syst. Sci., 28, 5107–5131,
+https://doi.org/10.5194/hess-28-5107-2024, <https://doi.org/10.5194/hess-28-5107-2024>.
+
 Schaller, N., Sillmann, J., Müller, M., Haarsma, R., Hazeleger, W., Hegdahl, T.J., Kelder, T.,
 van den Oord, G., Weerts, A., Whan, K., 2020. The role of spatial and temporal model resolution
 in a flood event storyline approach in western Norway. Weather and Climate Extremes, doi:

docs/model_docs/parameters_vertical.qmd

@@ -40,7 +40,7 @@ the `layered` profile, input parameter `kv` is used, and for the `layered_expone
 | **`whc`**  | water holding capacity as fraction of current snow pack  | - | 0.1  |
 | **`w_soil`**  | soil temperature smooth factor  | - | 0.1125  |
 | **`cf_soil`**  | controls soil infiltration reduction factor when soil is frozen  | - | 0.038  |
-| **`g_tt`**  | threshold temperature for glacier melt  | ᵒC | 0.0  |
+| **`g_ttm`**  | threshold temperature for glacier melt  | ᵒC | 0.0  |
 | **`g_cfmax`**  | Degree-day factor for glacier  | mm ᵒC$^{-1}$ Δt$^{-1}$| 3.0 mm ᵒC$^{-1}$ day$^{-1}$ |
 | **`g_sifrac`**  | fraction of the snowpack on top of the glacier converted into ice  | Δt$^{-1}$ | 0.001 day$^{-1}$ |
 | **`glacierfrac`**  | fraction covered by a glacier | - | 0.0  |

docs/model_docs/vertical/shared_processes.qmd

@@ -17,18 +17,18 @@ parameter `tti` defines how precipitation can occur partly as rain or snowfall.
 # https://gist.github.com/JoostBuitink/21dd32e71fd1360117fcd1c532c4fd9d#file-snowfall_fig-py # hide
 -->
 
-If precipitation occurs as snowfall, it is added to the dry snow component
-within the snow pack. Otherwise it ends up in the free water reservoir, which represents the
-liquid water content of the snow pack. Between the two components of the snow pack,
-interactions take place, either through snow melt (if temperatures are above a threshold `tt`)
-or through snow refreezing (if temperatures are below threshold `tt`).
+If precipitation occurs as snowfall, it is added to the dry snow component within the snow
+pack. Otherwise it ends up in the free water reservoir, which represents the liquid water
+content of the snow pack. Between the two components of the snow pack, interactions take
+place, either through snow melt (if temperatures are above a threshold `ttm`) or through
+snow refreezing (if temperatures are below threshold `ttm`).
 
 The respective rates of snow melt and refreezing are:
 
 $$
 \begin{align*}
-    Q_m &=& \subtext{\mathrm{cf}}{max}(T_a−\mathrm{tt})\, &&T_a > \mathrm{tt} \\
-    Q_r &=& \subtext{\mathrm{cf}}{max} \, \mathrm{cf}_r(\mathrm{tt}−T_a) &&T_a < \mathrm{tt}
+    Q_m &=& \subtext{\mathrm{cf}}{max}(T_a−\mathrm{ttm})\, &&T_a > \mathrm{ttm} \\
+    Q_r &=& \subtext{\mathrm{cf}}{max} \, \mathrm{cf}_r(\mathrm{ttm}−T_a) &&T_a < \mathrm{ttm}
 \end{align*}
 $$
 
@@ -61,34 +61,34 @@ into firn/ice (using the HBV-light model) and glacier melt (using a temperature
 model).
 
 The definition of glacier boundaries and initial volume is defined by two parameters. The
-parameter `glacierfrac` gives the fraction of each grid cell covered by a glacier as a number
-between zero and one. The state parameter `glacierstore` gives the amount of water (in mm w.e.)
-within the glaciers at each grid cell. Because the glacier store (`glacierstore`) cannot be
-initialized by running thFe model for a couple of years, a default initial state should be
-supplied by adding this parameter to the input static file. The required glacier data can be
-prepared from available glacier datasets.
+parameter `glacierfrac` gives the fraction of each grid cell covered by a glacier as a
+number between zero and one. The state parameter `glacierstore` gives the amount of water
+(in mm w.e.) within the glaciers at each grid cell. Because the glacier store
+(`glacierstore`) cannot be initialized by running the model for a couple of years, a default
+initial state should be supplied by adding this parameter to the input static file. The
+required glacier data can be prepared from available glacier datasets.
 
 First, a fixed fraction of the snowpack on top of the glacier is converted into ice for each
 timestep and added to the `glacierstore` using the HBV-light model (Seibert et al., 2018). This
 fraction `g_sifrac` typically ranges from $0.001$ to $0.006$.
 
 Then, when the snowpack on top of the glacier is almost all melted (snow cover $<
 \SI{10}{mm}$), glacier melt is enabled and estimated with a degree-day model. If the air
-temperature, $T_a$, is below a certain threshold `g_tt` ($\SIb{}{\degree C}$) precipitation
-occurs as snowfall, whereas it occurs as rainfall if $T_a ≥$ `g_tt`.
+temperature, $T_a$, is above a certain threshold `g_ttm` ($\SIb{}{\degree C}$) glacier melt
+occurs.
 
 With this the rate of glacier melt in mm is estimated as:
 
 $$
-Q_m = \subtext{g}{cfmax}(T_a − \subtext{g}{tt})\, ; \, T_a > \subtext{g}{tt}
+Q_m = \subtext{g}{cfmax}(T_a − \subtext{g}{ttm})\, ; \, T_a > \subtext{g}{ttm}
 $$
 
 where $Q_m$ is the rate of glacier melt and $\SIb{\subtext{g}{cfmax}}{mm (\degree
-C)^{-1}day^{-1}}$ is the melting factor. Parameter `g_tt` can be taken as equal to the snow
-`tt` parameter. Values of the melting factor `g_cfmax` normally varies from one glacier to
+C)^{-1}day^{-1}}$ is the melting factor. Parameter `g_ttm` can be taken as equal to the snow
+`ttm` parameter. Values of the melting factor `g_cfmax` normally varies from one glacier to
 another and some values are reported in the literature. `g_cfmax` can also be estimated by
-multiplying snow `cfmax` by a factor between 1 and 2, to take into account the higher albedo of
-ice compared to snow.
+multiplying snow `cfmax` by a factor between 1 and 2, to take into account the higher albedo
+of ice compared to snow.
 
 ## Rainfall interception
 Both the Gash and Rutter models are available to estimate rainfall interception by the

docs/user_guide/model_config.qmd

@@ -18,7 +18,7 @@ water_holding_capacity = "WHC"
 glacierstore = "wflow_glacierstore"
 glacierfrac = "wflow_glacierfrac"
 g_cfmax = "G_Cfmax"
-g_tt = "G_TT"
+g_ttm = "G_TT"
 g_sifrac = "G_SIfrac"
 
 [state.vertical]

server/test/client.jl

@@ -32,13 +32,13 @@ end
 
 @testset "model information functions" begin
     @test request((fn = "get_component_name",)) == Dict("component_name" => "sbm")
-    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 203)
-    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 203)
+    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 202)
+    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 202)
     to_check = [
         "vertical.soil.parameters.nlayers",
         "vertical.soil.parameters.theta_r",
-        "lateral.river.q",
-        "lateral.river.reservoir.outflow",
+        "lateral.river.variables.q",
+        "lateral.river.boundary_conditions.reservoir.variables.outflow",
     ]
     retrieved_vars = request((fn = "get_input_var_names",))["input_var_names"]
     @test all(x -> x in retrieved_vars, to_check)
@@ -49,12 +49,12 @@ end
 zi_size = 0
 vwc_1_size = 0
 @testset "variable information and get and set functions" begin
-    @test request((fn = "get_var_itemsize", name = "lateral.subsurface.ssf")) ==
+    @test request((fn = "get_var_itemsize", name = "lateral.subsurface.variables.ssf")) ==
           Dict("var_itemsize" => sizeof(Wflow.Float))
     @test request((fn = "get_var_type", name = "vertical.n"))["status"] == "ERROR"
     @test request((fn = "get_var_units", name = "vertical.soil.parameters.theta_s")) ==
           Dict("var_units" => "-")
-    @test request((fn = "get_var_location", name = "lateral.river.q")) ==
+    @test request((fn = "get_var_location", name = "lateral.river.variables.q")) ==
           Dict("var_location" => "node")
     zi_nbytes =
         request((fn = "get_var_nbytes", name = "vertical.soil.variables.zi"))["var_nbytes"]
@@ -68,19 +68,20 @@ vwc_1_size = 0
     vwc_1_itemsize =
         request((fn = "get_var_itemsize", name = "vertical.soil.variables.vwc[1]"))["var_itemsize"]
     vwc_1_size = Int(vwc_1_nbytes / vwc_1_itemsize)
-    @test request((fn = "get_var_grid", name = "lateral.river.h")) == Dict("var_grid" => 3)
+    @test request((fn = "get_var_grid", name = "lateral.river.variables.h")) ==
+          Dict("var_grid" => 3)
     msg = (fn = "get_value", name = "vertical.soil.variables.zi", dest = fill(0.0, zi_size))
     @test mean(request(msg)["value"]) ≈ 277.3620724821974
     msg = (fn = "get_value_ptr", name = "vertical.soil.parameters.theta_s")
     @test mean(request(msg)["value_ptr"]) ≈ 0.4409211971535584
     msg = (
         fn = "get_value_at_indices",
-        name = "lateral.river.q",
+        name = "lateral.river.variables.q",
         dest = [0.0, 0.0, 0.0],
         inds = [1, 5, 10],
     )
     @test request(msg)["value_at_indices"] ≈
-          [2.198747900215207f0, 2.6880427720508515f0, 3.4848783702629564f0]
+          [2.1007361866518766, 2.5702292750107687, 3.2904803551115727]
     msg =
         (fn = "set_value", name = "vertical.soil.variables.zi", src = fill(300.0, zi_size))
     @test request(msg) == Dict("status" => "OK")

server/test/sbm_config.toml

@@ -31,18 +31,18 @@ ustorelayerdepth = "ustorelayerdepth"
 snow_storage = "snow"
 snow_water = "snowwater"
 
-[state.lateral.river]
+[state.lateral.river.variables]
 h = "h_river"
 h_av = "h_av_river"
 q = "q_river"
 
-[state.lateral.river.reservoir]
+[state.lateral.river.boundary_conditions.reservoir.variables]
 volume = "volume_reservoir"
 
-[state.lateral.subsurface]
+[state.lateral.subsurface.variables]
 ssf = "ssf"
 
-[state.lateral.land]
+[state.lateral.land.variables]
 h = "h_land"
 h_av = "h_av_land"
 q = "q_land"
@@ -112,7 +112,7 @@ offset = 0.0
 
 [input.lateral.river]
 length = "wflow_riverlength"
-n = "N_River"
+mannings_n = "N_River"
 slope = "RiverSlope"
 width = "wflow_riverwidth"
 bankfull_elevation = "RiverZ"
@@ -132,7 +132,7 @@ targetminfrac = "ResTargetMinFrac"
 ksathorfrac = "KsatHorFrac"
 
 [input.lateral.land]
-n = "N"
+mannings_n = "N"
 slope = "Slope"
 
 [model]
@@ -161,17 +161,17 @@ ustorelayerdepth = "ustorelayerdepth"
 snow_storage = "snow"
 snow_water = "snowwater"
 
-[output.lateral.river]
+[output.lateral.river.variables]
 h = "h_river"
 q = "q_river"
 
-[output.lateral.river.reservoir]
+[output.lateral.river.boundary_conditions.reservoir.variables]
 volume = "volume_reservoir"
 
-[output.lateral.subsurface]
+[output.lateral.subsurface.variables]
 ssf = "ssf"
 
-[output.lateral.land]
+[output.lateral.land.variables]
 h = "h_land"
 q = "q_land"
 
@@ -181,7 +181,7 @@ path = "output_scalar_moselle.nc"
 [[netcdf.variable]]
 name = "Q"
 map = "gauges"
-parameter = "lateral.river.q"
+parameter = "lateral.river.variables.q"
 
 [[netcdf.variable]]
 coordinate.x = 6.255
@@ -202,13 +202,13 @@ path = "output_moselle.csv"
 
 [[csv.column]]
 header = "Q"
-parameter = "lateral.river.q"
+parameter = "lateral.river.variables.q"
 reducer = "maximum"
 
 [[csv.column]]
 header = "volume"
 index = 1
-parameter = "lateral.river.reservoir.volume"
+parameter = "lateral.river.boundary_conditions.reservoir.variables.volume"
 
 [[csv.column]]
 coordinate.x = 6.255
@@ -232,7 +232,7 @@ parameter = "vertical.atmospheric_forcing.temperature"
 [[csv.column]]
 header = "Q"
 map = "gauges"
-parameter = "lateral.river.q"
+parameter = "lateral.river.variables.q"
 
 [[csv.column]]
 header = "recharge"
@@ -255,8 +255,19 @@ components = [
   "vertical.snow.boundary_conditions",
   "vertical.snow.variables",
   "vertical.snow.parameters",
-  "lateral.subsurface",
-  "lateral.land",
-  "lateral.river",
-  "lateral.river.reservoir",
+  "lateral.subsurface.boundary_conditions",
+  "lateral.subsurface.variables",
+  "lateral.subsurface.parameters",
+  "lateral.subsurface.parameters.kh_profile",
+  "lateral.land.boundary_conditions",
+  "lateral.land.variables",
+  "lateral.land.variables.flow",
+  "lateral.land.parameters",
+  "lateral.river.variables",
+  "lateral.river.parameters",
+  "lateral.river.parameters.flow",
+  "lateral.river.boundary_conditions",
+  "lateral.river.boundary_conditions.reservoir.boundary_conditions",
+  "lateral.river.boundary_conditions.reservoir.parameters",
+  "lateral.river.boundary_conditions.reservoir.variables",
 ]

src/Wflow.jl

@@ -49,7 +49,7 @@ using Parameters: @with_kw
 using Polyester: @batch
 using ProgressLogging: @progress
 using StaticArrays: SVector, pushfirst, setindex
-using Statistics: mean, median, quantile!
+using Statistics: mean, median, quantile!, quantile
 using TerminalLoggers
 using TOML: TOML
 
@@ -141,8 +141,17 @@ Base.show(io::IO, m::Model) = print(io, "model of type ", typeof(m))
 
 include("forcing.jl")
 include("parameters.jl")
-include("flow.jl")
-include("horizontal_process.jl")
+include("groundwater/connectivity.jl")
+include("groundwater/aquifer.jl")
+include("groundwater/boundary_conditions.jl")
+include("routing/timestepping.jl")
+include("routing/subsurface.jl")
+include("routing/reservoir.jl")
+include("routing/lake.jl")
+include("routing/surface_kinwave.jl")
+include("routing/surface_local_inertial.jl")
+include("routing/surface_routing.jl")
+include("routing/routing_process.jl")
 include("vegetation/rainfall_interception.jl")
 include("vegetation/canopy.jl")
 include("snow/snow_process.jl")
@@ -154,7 +163,6 @@ include("soil/soil.jl")
 include("soil/soil_process.jl")
 include("sbm.jl")
 include("demand/water_demand.jl")
-include("reservoir_lake.jl")
 include("sbm_model.jl")
 include("sediment/erosion/erosion_process.jl")
 include("sediment/erosion/rainfall_erosion.jl")
@@ -170,9 +178,6 @@ include("sediment/sediment_transport/river_transport.jl")
 include("erosion.jl")
 include("sediment_flux.jl")
 include("sediment_model.jl")
-include("groundwater/connectivity.jl")
-include("groundwater/aquifer.jl")
-include("groundwater/boundary_conditions.jl")
 include("sbm_gwf_model.jl")
 include("utils.jl")
 include("bmi.jl")

src/bmi.jl

@@ -153,7 +153,7 @@ function BMI.get_var_grid(model::Model, name::String)
     s = split(name, "[")
     key = symbols(first(s))
     if exchange(param(model, key))
-        type = typeof(param(model, key[1:(end - 1)]))
+        type = typeof(param(model, key[1:2]))
         return if :reservoir in key
             0
         elseif :lake in key
@@ -162,9 +162,9 @@ function BMI.get_var_grid(model::Model, name::String)
             2
         elseif :river in key
             3
-        elseif type <: ShallowWaterLand && occursin("x", s[end])
+        elseif type <: LocalInertialOverlandFlow && occursin("x", s[end])
             4
-        elseif type <: ShallowWaterLand && occursin("y", s[end])
+        elseif type <: LocalInertialOverlandFlow && occursin("y", s[end])
             5
         else
             6
@@ -374,7 +374,7 @@ function BMI.get_grid_edge_nodes(model::Model, grid::Int, edge_nodes::Vector{Int
     m = div(n, 2)
     # inactive nodes (boundary/ghost points) are set at -999
     if grid == 3
-        nodes_at_edge = adjacent_nodes_at_link(network.river.graph)
+        nodes_at_edge = adjacent_nodes_at_edge(network.river.graph)
         nodes_at_edge.dst[nodes_at_edge.dst .== m + 1] .= -999
         edge_nodes[range(1, n; step = 2)] = nodes_at_edge.src
         edge_nodes[range(2, n; step = 2)] = nodes_at_edge.dst