Refactor sediment model (#463)

Fix issue #382.

---------

Co-authored-by: Willem van Verseveld <[email protected]>

build/create_binaries/download_test_data.jl

@@ -24,7 +24,7 @@ forcing_moselle_path = testdata(v"0.2.6", "forcing-moselle.nc", "forcing-moselle
 forcing_moselle_sed_path =
     testdata(v"0.2.3", "forcing-moselle-sed.nc", "forcing-moselle-sed.nc")
 staticmaps_moselle_sed_path =
-    testdata(v"0.2.3", "staticmaps-moselle-sed.nc", "staticmaps-moselle-sed.nc")
+    testdata(v"0.3.0", "staticmaps-moselle-sed.nc", "staticmaps-moselle-sed.nc")
 instates_moselle_sed_path =
     testdata(v"0.2", "instates-moselle-sed.nc", "instates-moselle-sed.nc")
 instates_moselle_path = testdata(v"0.2.6", "instates-moselle.nc", "instates-moselle.nc")
@@ -50,4 +50,5 @@ forcing_calendar_noleap_path =
 forcing_piave_path = testdata(v"0.2.9", "inmaps-era5-2010-piave.nc", "forcing-piave.nc")
 staticmaps_piave_path = testdata(v"0.2.9", "staticmaps-piave.nc", "staticmaps-piave.nc")
 instates_piave_path = testdata(v"0.2.9", "instates-piave.nc", "instates-piave.nc")
-instates_piave_gwf_path = testdata(v"0.2.9", "instates-piave-gwf.nc", "instates-piave-gwf.nc")
+instates_piave_gwf_path =
+    testdata(v"0.2.9", "instates-piave-gwf.nc", "instates-piave-gwf.nc")

docs/changelog.qmd

@@ -38,6 +38,11 @@ project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   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).
+- Refactor the `sediment` concept: similar to the refactor of `SBM` concept, the sediment
+  model has been split into three main structures. `SoilLoss` was split into rainfall_erosion,
+  overland_flow_erosion and (total) soil_erosion concepts. `OverlandFlowSediment` has been
+  split into transport_capacity, sediment_flux and to_river concepts. `RiverSediment` has been
+  split into transport_capacity, potential_erosion, sediment_flux and concentrations concepts.
 - 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

src/Wflow.jl

@@ -163,8 +163,20 @@ include("soil/soil.jl")
 include("soil/soil_process.jl")
 include("sbm.jl")
 include("demand/water_demand.jl")
-include("sediment.jl")
 include("sbm_model.jl")
+include("sediment/erosion/erosion_process.jl")
+include("sediment/erosion/rainfall_erosion.jl")
+include("sediment/erosion/overland_flow_erosion.jl")
+include("sediment/erosion/soil_erosion.jl")
+include("sediment/erosion/river_erosion.jl")
+include("sediment/sediment_transport/deposition.jl")
+include("sediment/sediment_transport/transport_capacity_process.jl")
+include("sediment/sediment_transport/transport_capacity.jl")
+include("sediment/sediment_transport/overland_flow_transport.jl")
+include("sediment/sediment_transport/land_to_river.jl")
+include("sediment/sediment_transport/river_transport.jl")
+include("erosion.jl")
+include("sediment_flux.jl")
 include("sediment_model.jl")
 include("sbm_gwf_model.jl")
 include("utils.jl")

src/erosion.jl

@@ -0,0 +1,81 @@
+"Soil loss model"
+@with_kw struct SoilLoss{RE, OFE, SE, T}
+    atmospheric_forcing::AtmosphericForcing{T}
+    hydrological_forcing::HydrologicalForcing{T}
+    geometry::LandGeometry{T}
+    rainfall_erosion::RE
+    overland_flow_erosion::OFE
+    soil_erosion::SE
+end
+
+"Initialize soil loss model"
+function SoilLoss(dataset, config, indices)
+    n = length(indices)
+
+    atmospheric_forcing = AtmosphericForcing(n)
+    hydrological_forcing = HydrologicalForcing(n)
+    geometry = LandGeometry(dataset, config, indices)
+
+    # Rainfall erosion
+    rainfallerosionmodel = get(config.model, "rainfall_erosion", "answers")::String
+    if rainfallerosionmodel == "answers"
+        rainfall_erosion_model = RainfallErosionAnswersModel(dataset, config, indices)
+    elseif rainfallerosionmodel == "eurosem"
+        rainfall_erosion_model = RainfallErosionEurosemModel(dataset, config, indices)
+    else
+        error("Unknown rainfall erosion model: $rainfallerosionmodel")
+    end
+
+    # Overland flow erosion
+    overlandflowerosionmodel = get(config.model, "overland_flow_erosion", "answers")::String
+    if overlandflowerosionmodel == "answers"
+        overland_flow_erosion_model =
+            OverlandFlowErosionAnswersModel(dataset, config, indices)
+    else
+        error("Unknown overland flow erosion model: $overlandflowerosionmodel")
+        # overland_flow_erosion_model = NoOverlandFlowErosionModel()
+    end
+
+    # Total soil erosion and particle differentiation
+    soil_erosion_model = SoilErosionModel(dataset, config, indices)
+
+    soil_loss = SoilLoss{
+        typeof(rainfall_erosion_model),
+        typeof(overland_flow_erosion_model),
+        typeof(soil_erosion_model),
+        Float,
+    }(;
+        atmospheric_forcing = atmospheric_forcing,
+        hydrological_forcing = hydrological_forcing,
+        geometry = geometry,
+        rainfall_erosion = rainfall_erosion_model,
+        overland_flow_erosion = overland_flow_erosion_model,
+        soil_erosion = soil_erosion_model,
+    )
+    return soil_loss
+end
+
+"Update soil loss model for a single timestep"
+function update!(model::SoilLoss, dt)
+    (;
+        atmospheric_forcing,
+        hydrological_forcing,
+        geometry,
+        rainfall_erosion,
+        overland_flow_erosion,
+        soil_erosion,
+    ) = model
+
+    #TODO add interception/canopygapfraction calculation here for eurosem
+    #need SBM refactor
+
+    # Rainfall erosion
+    update_boundary_conditions!(rainfall_erosion, atmospheric_forcing, hydrological_forcing)
+    update!(rainfall_erosion, geometry, dt)
+    # Overland flow erosion
+    update_boundary_conditions!(overland_flow_erosion, hydrological_forcing)
+    update!(overland_flow_erosion, geometry, dt)
+    # Total soil erosion and particle differentiation
+    update_boundary_conditions!(soil_erosion, rainfall_erosion, overland_flow_erosion)
+    update!(soil_erosion)
+end

src/forcing.jl

@@ -1,5 +1,5 @@
 "Struct to store atmospheric forcing variables"
-@get_units @with_kw struct AtmosphericForcing{T}
+@get_units @grid_loc @with_kw struct AtmosphericForcing{T}
     # Precipitation [mm Δt⁻¹]
     precipitation::Vector{T}
     # Potential reference evapotranspiration [mm Δt⁻¹]
@@ -15,9 +15,37 @@ function AtmosphericForcing(
     potential_evaporation::Vector{T} = fill(mv, n),
     temperature::Vector{T} = fill(mv, n),
 ) where {T}
-    return AtmosphericForcing{T}(;
-        precipitation,
-        potential_evaporation,
-        temperature,
+    return AtmosphericForcing{T}(; precipitation, potential_evaporation, temperature)
+end
+
+"Struct to store hydrological forcing variables"
+@get_units @grid_loc @with_kw struct HydrologicalForcing{T}
+    # Rainfall interception by the vegetation [mm]
+    interception::Vector{T} | "mm"
+    # Overland flow depth [m]
+    waterlevel_land::Vector{T} | "m"
+    # Overland flow discharge [m3 s-1]
+    q_land::Vector{T} | "m3 s-1"
+    # River depth [m]
+    waterlevel_river::Vector{T} | "m"
+    # River discharge [m3 s-1]
+    q_river::Vector{T} | "m3 s-1"
+end
+
+"Initialize hydrological forcing"
+function HydrologicalForcing(
+    n;
+    interception::Vector{T} = fill(mv, n),
+    waterlevel_land::Vector{T} = fill(mv, n),
+    q_land::Vector{T} = fill(mv, n),
+    waterlevel_river::Vector{T} = fill(mv, n),
+    q_river::Vector{T} = fill(mv, n),
+) where {T}
+    return HydrologicalForcing{T}(;
+        interception,
+        waterlevel_land,
+        q_land,
+        waterlevel_river,
+        q_river,
     )
 end

src/parameters.jl

@@ -101,4 +101,87 @@ function VegetationParameters(dataset, config, indices)
         )
     end
     return vegetation_parameter_set
+end
+
+"Struct to store land geometry parameters"
+@get_units @grid_loc @with_kw struct LandGeometry{T}
+    # cell area [m^2]
+    area::Vector{T} | "m^2"
+    # drain width [m]
+    width::Vector{T} | "m"
+    # drain slope
+    slope::Vector{T} | "-"
+end
+
+"Initialize land geometry parameters"
+function LandGeometry(nc, config, inds)
+    # read x, y coordinates and calculate cell length [m]
+    y_nc = read_y_axis(nc)
+    x_nc = read_x_axis(nc)
+    y = permutedims(repeat(y_nc; outer = (1, length(x_nc))))[inds]
+    cellength = abs(mean(diff(x_nc)))
+
+    sizeinmetres = get(config.model, "sizeinmetres", false)::Bool
+    xl, yl = cell_lengths(y, cellength, sizeinmetres)
+    area = xl .* yl
+    ldd = ncread(nc, config, "ldd"; optional = false, sel = inds, allow_missing = true)
+    drain_width = map(get_flow_width, ldd, xl, yl)
+    landslope = ncread(
+        nc,
+        config,
+        "vertical.land_parameter_set.slope";
+        optional = false,
+        sel = inds,
+        type = Float,
+    )
+    clamp!(landslope, 0.00001, Inf)
+
+    land_parameter_set =
+        LandGeometry{Float}(; area = area, width = drain_width, slope = landslope)
+    return land_parameter_set
+end
+
+"Struct to store river geometry parameters"
+@get_units @grid_loc @with_kw struct RiverGeometry{T}
+    # river width [m]
+    width::Vector{T} | "m"
+    # river length
+    length::Vector{T} | "m"
+    # slope
+    slope::Vector{T} | "-"
+end
+
+"Initialize river geometry parameters"
+function RiverGeometry(nc, config, inds)
+    riverwidth = ncread(
+        nc,
+        config,
+        "lateral.river_parameter_set.width";
+        optional = false,
+        sel = inds,
+        type = Float,
+    )
+    riverlength = ncread(
+        nc,
+        config,
+        "lateral.river_parameter_set.length";
+        optional = false,
+        sel = inds,
+        type = Float,
+    )
+    riverslope = ncread(
+        nc,
+        config,
+        "lateral.river_parameter_set.slope";
+        optional = false,
+        sel = inds,
+        type = Float,
+    )
+    minimum(riverlength) > 0 || error("river length must be positive on river cells")
+    minimum(riverwidth) > 0 || error("river width must be positive on river cells")
+    clamp!(riverslope, 0.00001, Inf)
+
+    river_parameter_set =
+        RiverGeometry{Float}(; width = riverwidth, length = riverlength, slope = riverslope)
+    return river_parameter_set
 end