Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor sediment model #463

Merged
merged 45 commits into from
Dec 9, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
45 commits
Select commit Hold shift + click to select a range
de91ed9
Remove vertical concepts `HBV` and `FLEXTopo` (#433)
vers-w Jul 1, 2024
5107824
Cleanup metadata macros (#434)
vers-w Jul 2, 2024
d21c930
Refactor/style guide (#437)
CFBaptista Jul 15, 2024
5ab02bf
Stop using local JULIAUP_DEPOT_PATH (#438) (#440)
vers-w Jul 18, 2024
9819e21
Merge branch 'master' into v1
vers-w Sep 4, 2024
5e87215
first commit refactor soil erosion
hboisgon Sep 9, 2024
a98b2a4
small updates and bugfixes until test comparison
hboisgon Sep 17, 2024
0288341
Working tests for soil erosion
hboisgon Sep 17, 2024
38f90e4
re-enable all tests and remove non-used config
hboisgon Sep 17, 2024
3847756
refactor sediment lateral land
hboisgon Sep 19, 2024
4f0383e
working refactor of lateral land sediment
hboisgon Sep 20, 2024
b241304
draft refactor river sediment
hboisgon Sep 23, 2024
903508c
working tests and all results comparable to master
hboisgon Sep 25, 2024
e932293
renaming of functions
hboisgon Sep 25, 2024
55b48e9
Remove vertical concepts `HBV` and `FLEXTopo` (#433)
vers-w Jul 1, 2024
e1adf77
Cleanup metadata macros (#434)
vers-w Jul 2, 2024
4181a6d
Refactor/style guide (#437)
CFBaptista Jul 15, 2024
992671a
first commit refactor soil erosion
hboisgon Sep 9, 2024
39d938c
small updates and bugfixes until test comparison
hboisgon Sep 17, 2024
4fbb8b5
Working tests for soil erosion
hboisgon Sep 17, 2024
3db93fb
re-enable all tests and remove non-used config
hboisgon Sep 17, 2024
f4f7124
refactor sediment lateral land
hboisgon Sep 19, 2024
1014297
working refactor of lateral land sediment
hboisgon Sep 20, 2024
ed03716
draft refactor river sediment
hboisgon Sep 23, 2024
15ff556
working tests and all results comparable to master
hboisgon Sep 25, 2024
4db4d48
renaming of functions
hboisgon Sep 25, 2024
4cd7a23
Merge branch 'refactor/sediment_concept' of https://github.com/Deltar…
hboisgon Nov 20, 2024
f5191c2
undo wrong changes from rebase
hboisgon Nov 20, 2024
c7fa59b
Update BMI Sediment
vers-w Nov 26, 2024
a3ffec7
Remove unpack
vers-w Nov 26, 2024
c15cf6a
Remove duplicate `get_sediment_states()`
vers-w Nov 26, 2024
c53fa36
Update version testdata `staticmaps-moselle-sed.nc`
vers-w Nov 26, 2024
b2ff488
Rename couple of initialize functions Sediment
vers-w Nov 26, 2024
6b62c34
Add BMI `@grid_loc` macro to erosion files
vers-w Nov 26, 2024
03c57b0
move geometry to parameters
hboisgon Dec 4, 2024
faaf49b
move sediment files to a sediment folder
hboisgon Dec 4, 2024
70e3516
split atmospheric and hydrological forcing
hboisgon Dec 4, 2024
ab48727
rename nc, inds to dataset, indices
hboisgon Dec 4, 2024
3b99f49
rename ts to dt
hboisgon Dec 4, 2024
ff673c3
add docstrings
hboisgon Dec 5, 2024
36c155e
improve test coverage
hboisgon Dec 5, 2024
22c2f45
undo skipping tests
hboisgon Dec 5, 2024
67ffcb2
Merge branch 'master' into refactor/sediment_concept
hboisgon Dec 5, 2024
2246a3d
bugfixes after merge + changelog
hboisgon Dec 5, 2024
af9dd7a
use geometry again
hboisgon Dec 6, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions build/create_binaries/download_test_data.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand All @@ -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")
5 changes: 5 additions & 0 deletions docs/changelog.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
14 changes: 13 additions & 1 deletion src/Wflow.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
81 changes: 81 additions & 0 deletions src/erosion.jl
Original file line number Diff line number Diff line change
@@ -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
38 changes: 33 additions & 5 deletions src/forcing.jl
Original file line number Diff line number Diff line change
@@ -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⁻¹]
Expand All @@ -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
83 changes: 83 additions & 0 deletions src/parameters.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading
Loading