Refactor SBM concept (#487)

* Remove vertical concepts `HBV` and `FLEXTopo` (#433)

* Remove vertical concept `FLEXTopo`

* Remove vertical concept `HBV`

* Update download test data for build

* Cleanup docs
Removed `HBV` and `FLEXTopo` concepts.

* Removed code related to `FLEXTopo`
Use of extra dim `classes`.

* Remove and change river and land `inwater` functions
As a result of removing `HBV` and `FLEXTopo` concepts.

* Update changelog

* Fix typo in docs

Co-authored-by: JoostBuitink <[email protected]>

---------

Co-authored-by: JoostBuitink <[email protected]>

* Cleanup metadata macros (#434)

* Cleanup metadata structs
Remove `exchange`, `grid_type` and `grid_location` from metadata structs.  `grid_type` is not required, it is already implemented as part of BMI. `exhange` and `grid_location` are now implemented without the use of the FieldMetadata package.

* Update changelog

* Fix `exchange` function for `ShallowWaterLand`

* Add tests

* Move `exchange` and `grid_location` functions
To file bmi.jl as these functions are quite specific to BMI functionality.

* Add `v1` branch to CI

* Refactor/style guide (#437)

* Add configuration file for Julia Formatter

* Add format on save in vscode settings file

* Remove .vscode from gitignore

* Format all Julia files

* Sync wflow's style guide with ribasim's style guide

* Update formatting

* Update .vscode/settings.json

* Remove outdated comments

* Add workaround for VS Code bugs

* WIP: refactor `SBM` interception

* WIP: add `AtmosphericForcing` to store meteo input

* WIP: remove unused type parameter `A` from `SBM`

* WIP: add `RutterInterceptionModel`

* WIP: refactor `SBM` snow

* Stop using local JULIAUP_DEPOT_PATH (#438) (#440)

* Stop using local JULIAUP_DEPOT_PATH

* Avoid error when override is already set

Co-authored-by: Martijn Visser <[email protected]>

* WIP: refactor `SBM` glacier

* WIP: introduce types `NoSnowModel` and `NoGlacierModel`

* WIP:  `model` in `update` functions and add wrapper methods
In `update` functions use `model` for interception, snow and glacier models. Add wrapper methods for snow runoff and glacier melt.

* WIP: start refactor soil model of `SBM` (initialization)

* Update of `SBM` soil model (until recharge)

* Update run of refactored sbm concept (sbm_model)
This is with kinematic wave routing. Renamed the model for `SBM` from `soil` to `bucket` as it is more related to this type of model (Simple Bucket Model) than a soil model that for example solves the Richards equation for the unsaturated zone and upper part of the saturated zone.

* Delete `soil` folder

* Add `SBM` functions to seperate file (bucket_process.jl)
Also removed the vertical_process.jl file.

* Add water demand functionality

* Rename `vegetation_parameters` to `vegetation_parameter_set`

* Fix test Wflow ZMQ Server

* Refactor update of `LandHydrologySBM`
Add water demand and allocation, and snow transport computations.

* Further refactor of update `LandHydrologySBM`
Add separate functions to update the boundary conditions of the snow and `SBM` models.

* Rename `bucket` to `soil` and add `SurfaceRunoff` model
In addition:
- add a parameter set for LandParameters (to share accross model components).
- make use of type NamedTuple (for external models) in the update of models and boundary conditions, these external models provide input to the models and boundary conditions.
- move shared parameters (vegetation and land) to file Parameters.jl.

* Include `atmospheric_forcing` as a separate argument
In update of models and boundaries of `LandHydrologySBM`

* Remove land_parameter_set
It only contains two parameters that are not shared that much across different models of `LandHydrologySBM`.

* Rename outer constructors and couple of structs
Use the struct name for outer constructor functions.

* Add parameter `soil_fraction` to `SbmSoilParameters`

* Move function for `soil_fraction` to soil.jl

* Rename model  type `SurfaceRunoff` to `OpenWaterRunoff`

* Refactor water demand
For more efficient and cleaner code: 1) Support dispatching on types that represent non-existent Demand (e.g. NoDemand) or Allocation (e.g. NoAllocationLand) models is supported. 2) Add more wrapper functions to access variables of water demand and allocation structs.

* Move runoff,jl to separate folder

* Refactor update `SbmSoilModel`
Add functions for different processes that are part of the `SbmSoilModel` to improve code readability.

* Use types for hydraulic conductivity profiles
And reduce allocations.

* Rename `surface_water_flux` to `water_flux_surface`

* Small updates

* Boundary conditions `SbmSoilModel`
Change computation of `water_flux_surface` and `potential_soilevaporation` to represent fluxes for actual soil boundary by moving paddy evaporation and infiltration to the `update_boundary_conditions!` function.

* FIx some non-concrete types

* Delete log.txt file

* Refactor water demand
Split structs into variables and parameters and moved demand variables of `AllocationLand` to `Demand`.

* Add docstrings to snow model

* Add docstrings glacier model

* Add docstrings canopy and rename `interception_flux` variable

* Add docstrings open water runoff model

* Add docstrings SBM soil model

* Add docstrings water demand
And group structs and associated functions.

* Add comment about metadata and BMI (Wflow.jl)

* Add docstrings to vegetation parameters

* Add docstrings atmospheric forcing

* Cleanup metadata `get_units` and `grid_loc`
Only specify for variables that are exposed by BMI. Additionally added docstrings for `LandHydrologySBM` and simplified initialization of `LandHydrologySBM`.

* Update docstrings and comments

* Fix ZMQ Server tests

* Ksat_profile tests from PR #474
So we can close PR #474 as the fix is already part if this branch.

* Update changelog

* Update docstrings

* Add variable `snow_melt` to `SnowVariables`

* Update changelog

* Format Wflow.jl

* Suggestions from code review for water_demand

Co-authored-by: Bart de Koning <[email protected]>

* Address remaining review comments water_demand

* Remaining suggestions from code review

Co-authored-by: Bart de Koning <[email protected]>

* Address review comments
Simplify external constructors.

* Address review comment
Rewrite without explicit loop.

* Address review comment and bug fix
-Use consistently `!` for functions that are mutating, and return always  `nothing` for these functions.
- Fix missing `canonicalize` from `Dates`

* Fix typo

Co-authored-by: JoostBuitink <[email protected]>

* Address remaining comments 2nd reviewer

* Fix wflow server tests

* Revert back to using `isnothing`
As of Julia 1.7 there is no difference and `isnothing` is also used in other code sections.

---------

Co-authored-by: JoostBuitink <[email protected]>
Co-authored-by: Carlos Fernando Baptista <[email protected]>
Co-authored-by: Martijn Visser <[email protected]>
Co-authored-by: Bart de Koning <[email protected]>

.JuliaFormatter.toml

@@ -0,0 +1,12 @@
+# Options for the JuliaFormatter auto syntax formatting tool.
+# https://domluna.github.io/JuliaFormatter.jl/stable/
+# https://docs.sciml.ai/SciMLStyle/stable/
+
+# Based on the default style we do pick these non-default options from SciML style:
+whitespace_ops_in_indices = true
+remove_extra_newlines = true
+always_for_in = true
+whitespace_typedefs = true
+
+# And add other options we like:
+separate_kwargs_with_semicolon = true

.github/workflows/CI.yml

@@ -4,6 +4,7 @@ on:
   push:
     branches:
       - master
+      - v1
     tags: '*'
 jobs:
   test:

.github/workflows/CIWflowServer.yml

@@ -4,6 +4,7 @@ on:
   push:
     branches:
       - master
+      - v1
     tags: '*'
 jobs:
   test:

.gitignore

@@ -1,9 +1,6 @@
 # General
 .DS_Store
 
-# VS Code stuff
-.vscode
-
 # Packaging stuff
 *.jl.*.cov
 *.jl.cov

.vscode/settings.json

@@ -0,0 +1,9 @@
+{
+    "[julia]": {
+        "editor.formatOnSave": true,
+    },
+    "julia.lint.disabledDirs": [
+        ".pixi"
+    ],
+    "julia.lint.run": true,
+}

build/create_binaries/download_test_data.jl

@@ -19,12 +19,8 @@ end
 
 staticmaps_rhine_path = testdata(v"0.1", "staticmaps.nc", "staticmaps-rhine.nc")
 staticmaps_moselle_path =
-    testdata(v"0.2.8", "staticmaps-moselle.nc", "staticmaps-moselle.nc")
-staticmaps_lahn_path = testdata(v"0.2.1", "staticmaps-lahn.nc", "staticmaps-lahn.nc")
-staticmaps_meuse_path =
-    testdata(v"0.2.8", "staticmaps_flex_meuse.nc", "staticmaps_flex_meuse.nc")
+    testdata(v"0.2.9", "staticmaps-moselle.nc", "staticmaps-moselle.nc")
 forcing_moselle_path = testdata(v"0.2.6", "forcing-moselle.nc", "forcing-moselle.nc")
-forcing_lahn_path = testdata(v"0.2", "forcing-lahn.nc", "forcing-lahn.nc")
 forcing_moselle_sed_path =
     testdata(v"0.2.3", "forcing-moselle-sed.nc", "forcing-moselle-sed.nc")
 staticmaps_moselle_sed_path =
@@ -42,7 +38,6 @@ forcing_sbm_gw_path = testdata(
     "forcing-sbm-groundwater-part2.nc",
     "forcing-sbm-groundwater-part2.nc",
 )
-forcing_meuse_path = testdata(v"0.2.8", "forcing_meuse.nc", "forcing_meuse.nc")
 staticmaps_sbm_gw_path =
     testdata(v"0.2.3", "staticmaps-sbm-groundwater.nc", "staticmaps-sbm-groundwater.nc")
 instates_sbm_gw_path =

build/wflow_cli/test/runtests.jl

@@ -2,7 +2,8 @@ using Test
 using Wflow
 
 extension = Sys.iswindows() ? ".exe" : ""
-wflow_exe = normpath(@__DIR__, "../../create_binaries/wflow_bundle/bin/wflow_cli" * extension)
+wflow_exe =
+    normpath(@__DIR__, "../../create_binaries/wflow_bundle/bin/wflow_cli" * extension)
 
 # this assumes that the Wflow tests have already been run, so the data has been downloaded
 testdir = abspath(dirname(pathof(Wflow)), "..", "test")
@@ -14,15 +15,15 @@ outputdir = joinpath(datadir, "output")
 
 @testset "no_config" begin
     # Clean output directory
-    rm(outputdir; force=true, recursive=true)
+    rm(outputdir; force = true, recursive = true)
     @test_throws ProcessFailedException run(`$wflow_exe`)
     # Check if no files are being created
     @test !(isdir(outputdir))
 end
 
 @testset "wflow_sbm" begin
     # Clean directory
-    rm(outputdir; force=true, recursive=true)
+    rm(outputdir; force = true, recursive = true)
     # Run cli with the toml
     toml = normpath(testdir, "sbm_config.toml")
     run(`$wflow_exe $toml`)
@@ -43,7 +44,7 @@ end
 
 @testset "wflow_sbm-gwf" begin
     # Clean directory
-    rm(outputdir; force=true, recursive=true)
+    rm(outputdir; force = true, recursive = true)
     # Run cli with the toml
     toml = normpath(testdir, "sbm_gwf_config.toml")
     run(`$wflow_exe $toml`)
@@ -63,7 +64,7 @@ end
 
 @testset "wflow_sediment" begin
     # Clean directory
-    rm(outputdir; force=true, recursive=true)
+    rm(outputdir; force = true, recursive = true)
     # Run cli with the toml
     toml = normpath(testdir, "sediment_config.toml")
     run(`$wflow_exe $toml`)
@@ -82,11 +83,12 @@ end
 end
 
 @testset "wflow_sbm_timing" begin
-
     toml = normpath(testdir, "sbm_config.toml")
-    time_sbm_1thread = @elapsed run(Cmd(`$wflow_exe $toml`, env=("JULIA_NUM_THREADS" => "1",)))
+    time_sbm_1thread =
+        @elapsed run(Cmd(`$wflow_exe $toml`; env = ("JULIA_NUM_THREADS" => "1",)))
 
-    time_sbm_4thread = @elapsed run(Cmd(`$wflow_exe $toml`, env=("JULIA_NUM_THREADS" => "4", )))
+    time_sbm_4thread =
+        @elapsed run(Cmd(`$wflow_exe $toml`; env = ("JULIA_NUM_THREADS" => "4",)))
 
     # Test if run with more threads is indeed faster
     @test time_sbm_4thread < time_sbm_1thread

docs/make.jl

@@ -23,8 +23,6 @@ pages = [
         "model_docs/model_configurations.md",
         "Vertical concepts" => [
             "model_docs/vertical/sbm.md",
-            "model_docs/vertical/hbv.md",
-            "model_docs/vertical/flextopo.md",
             "model_docs/vertical/sediment.md",
             "model_docs/shared_concepts.md",
         ],

docs/src/changelog.md

@@ -5,6 +5,32 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## Unreleased
+
+### Fixed
+- Initialization of `LateralSSF` variables `ssf` and `ssfmax` with vertical hydraulic
+  conductivity profile `exponential_constant`. Removed parameter `khfrac` from the
+  computation, as it is already part of parameter `kh_0`.
+
+### Changed
+- Removed vertical concepts `HBV` and `FLEXTopo`.
+- Removed metadata macros `exchange` and `grid_type`. The macro `grid_type` is not required
+  because this functionality is already part of `BMI`. The macro `exchange` is replaced by a
+  function used by `BMI`. Remaining metadata macros `get_units` and `grid_loc` are only used
+  by `BMI`.
+- Refactor the vertical `SBM` concept: divide the long struct `SBM` into different model
+  components for interception, snow, glacier, (open water) runoff, soil, water demand and
+  allocation stored in the struct `LandHydrologySBM`. Additionally, the atmospheric forcing
+  and a shared vegetation parameterset are stored as separate fields in struct
+  `LandHydrologySBM` (with soil model `SbmSoilModel`). The model component structs have
+  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. 
+
+### Added
+- Support direct output of snow and glacier melt, and add computation of snow water
+  equivalent (SWE).
+
 ## v0.8.1 - 2024-08-27
 
 ### Fixed

docs/src/model_docs/model_configurations.md

@@ -131,128 +131,6 @@ lateral.land => struct ShallowWaterLand{T}
 The local inertial approach is described in more detail in the section [Local inertial
 model](@ref local_inertial).
 
-## [wflow\_hbv](@id config_hbv)
-The Hydrologiska Byrans Vattenbalansavdelning (HBV) model was introduced back in 1972 by the
-Swedish Meteological and Hydrological Institute (SMHI).  The HBV model is mainly used for
-runoff simulation and hydrological forecasting. The model is particularly useful for
-catchments where snow fall and snow melt are dominant factors, but application of the model
-is by no means restricted to these type of catchments.
-
-The model is based on the HBV-96 model. However, the hydrological routing represented in HBV
-by a triangular function controlled by the MAXBAS parameter has been removed. Instead, the
-kinematic wave function is used to route the water downstream. All runoff that is generated
-in a cell in one of the HBV reservoirs is added to the kinematic wave reservoir at the end
-of a timestep. There is no connection between the different HBV cells within the model.
-
-A catchment is divided into a number of grid cells. For each of the cells individually,
-daily runoff is computed through application of the HBV-96 of the HBV model. The use of the
-grid cells offers the possibility to turn the HBV modelling concept, which is originally
-lumped, into a distributed model.
-
-![wflow_hbv model](../images/hbv96.png)
-
-The figure above shows a schematic view of hydrological response simulation with the
-HBV-modelling concept. The land-phase of the hydrological cycle is represented by three
-different components: a snow routine, a soil routine and a runoff response routine. Each
-component is discussed in more detail below.
-
-The vertical HBV concept is described in section [HBV vertical concept](@ref vert_hbv). The
-routing for river and overland flow is described in the section [Kinematic wave](@ref
-kin_wave).
-
-Below the mapping for wflow\_hbv (type `hbv`) to the vertical HBV concept (instance of
-`struct HBV`) and the different lateral concepts. For an explanation about the type
-parameters between curly braces after the `struct` name see the section on model parameters.
-
-```julia
-vertical => struct HBV{T}
-lateral.subsurface => struct LateralSSF{T}
-lateral.land => struct SurfaceFlow{T,R,L}
-lateral.river => struct SurfaceFlow{T,R,L}
-lateral.river.lake => struct NaturalLake{T} # optional
-lateral.river.reservoir => struct SimpleReservoir{T} # optional
-```
-
-## [wflow\_flextopo](@id config_flextopo)
-The FLEXTopo model is a process-based model, which consists of different parallel classes
-connected through their groundwater storage. These classes are usually delineated from
-topographical data to represent the variability in hydrological processes across
-user-defined Hydrological Response Units (HRU). The main assumption underlying the concept,
-which was first introduced by Savenije (2010), is that different parts of the landscape
-fulfill different tasks in runoff generation and, hence, can be represented by different
-model structures. The strength of the concept is that the definition of classes and
-associated model structures is modular and flexible and not fixed to a predefined model
-structure. The flexible approach allows to develop process-based models for different
-topographic, climatic, geologic and land use conditions, making use of the available data
-and expert knowledge.
-
-The kinematic wave function is used to route the water downstream. In a similar way as for
-HBV, all runoff that is generated in a cell in one of the FLEXTopo storages is added to the
-kinematic wave reservoir at the end of a timestep. There is no connection between the
-different vertical FLEXTopo cells within the model. The FLEXTopo model is implemented in a
-fully distributed way in the wflow Julia framework.
-
-In wflow\_flextopo, the user is free to determine the number of classes and which model
-components to include or exclude for each class, this is done in the TOML file. Currently,
-for each storage, it is possible to bypass the storage and pass on the fluxes to the next
-model component. Interested users can contribute to the code by adding other
-conceptualizations for each storage components.
-
-```julia
-[model]
-type = "flextopo"
-classes = ["h", "p", "w"]     #user can set the number and name of each class.
-
-# for each component which is class specific, the user can select which conceptualization
-# to apply for each class as defined above in classes = ["h", "p", "w"]
-select_snow = ["common_snow_hbv"]
-# available options are ["common_snow_hbv", "common_snow_no_storage"]
-select_interception = ["interception_overflow", "interception_overflow", "interception_overflow"]
-# available options are ["interception_overflow", "interception_no_storage"]
-select_hortonponding = ["hortonponding_no_storage", "hortonponding_no_storage", "hortonponding_no_storage"]
-# available options are ["hortonponding", "hortonponding_no_storage"]
-select_hortonrunoff = ["hortonrunoff_no_storage", "hortonrunoff_no_storage", "hortonrunoff_no_storage"]
-# available options are ["hortonrunoff", "hortonrunoff_no_storage"]
-select_rootzone = ["rootzone_storage", "rootzone_storage", "rootzone_storage"]
-# available options are ["rootzone_storage", "rootzone_no_storage"]
-select_fast = ["fast_storage", "fast_storage", "fast_storage"]
-# available options are ["fast_storage", "fast_no_storage"]
-select_slow = ["common_slow_storage"]
-# available options are ["common_slow_storage", "slow_no_storage"]
-```
-
-A schematic representation of the most complete model structure including all storage
-components, as currently implemented in the code, is shown in the Figure below. When setting
-up the model with multiple classes, model structures can be adapted by bypassing storages or
-turning parameter values on or off (e.g.: percolation or capillary rise, non-linear versus
-linear outflow of the fast runoff etc.), an example of a three class model is shown in
-[FLEXTopo vertical concept](@ref vert_flextopo).
-
-![flextopo_julia_1class.png](../images/flextopo_julia_1class.png)
-*Schematic representation of the FLEXTopo model for a single class model including all
-storages and fluxes. Main parameters are denoted in red.*
-
-In the staticmaps, the user needs to provide maps of the fraction of each class within each
-cell, as shown below with `hrufrac`. For each model parameter which is class specific, an
-extra dimension `classes` is required in the staticmaps netcdf. For an example model, see
-[FLEXTopo example model](@ref wflow_flextopo_data).
-
-```julia
-[input.vertical]
-hrufrac = "hrufrac_lu"
-```
-
-Parameter multiplication of model parameters which are defined for several classes is
-possible through the TOML file:
-
-```julia
-[input.vertical.kf]
-netcdf.variable.name = "kf"
-scale = [1.0, 3.0, 4.0]
-offset = [0.0, 0.0, 0.0]
-class = ["h", "p", "w"]
-```
-
 ## [wflow\_sediment](@id config_sediment)
 The processes and fate of many particles and pollutants impacting water quality at the
 catchment level are intricately linked to the processes governing sediment dynamics. Both
@@ -292,8 +170,7 @@ required dynamic inputs to run wflow\_sediment are:
 -  River water level in the kinematic wave,
 -  Rainfall interception by the vegetation.
 
-These inputs can be obtained from other wflow models such as wflow\_sbm, wflow\_hbv or from
-other sources.
+These inputs can be obtained from wflow\_sbm or from other sources.
 
 Model outputs can be saved for both the inland and the instream part of the model. Some
 examples are listed below.
@@ -331,7 +208,4 @@ Bedconc = "Bedconc"
 ## References
 + Köhler, L., Mulligan, M., Schellekens, J., Schmid, S., Tobón, C., 2006, Hydrological
   impacts of converting tropical montane cloud forest to pasture, with initial reference to
-  northern Costa Rica. Final Technical Report DFID‐FRP Project No. R799.
-
-+ Savenije, H. H. G. (2010). HESS opinions “topography driven conceptual modelling (FLEX-Topo).”
-  Hydrology and Earth System Sciences, 14(12), 2681–2692. https://doi.org/10.5194/hess-14-2681-2010
+  northern Costa Rica. Final Technical Report DFID‐FRP Project No. R799.