From 490c606d8d52eee7f58631c64ce7cc96b768dc6d Mon Sep 17 00:00:00 2001 From: Willem van Verseveld Date: Thu, 22 May 2025 15:03:48 +0200 Subject: [PATCH 1/5] Fix init state local inertial land routing --- Wflow/src/sbm_model.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/Wflow/src/sbm_model.jl b/Wflow/src/sbm_model.jl index 5cbb422b0..e310eeb7a 100644 --- a/Wflow/src/sbm_model.jl +++ b/Wflow/src/sbm_model.jl @@ -168,11 +168,11 @@ function set_states!(model::AbstractModel{<:Union{SbmModel, SbmGwfModel}}) end land_v.storage .= land_v.h .* surface_flow_width .* flow_length elseif land_routing == "local-inertial" - (; river, x_length, y_length) = domain.land.parameters + (; river_location, x_length, y_length) = domain.land.parameters (; flow_width, flow_length) = domain.river.parameters - (; bankfull_storage) = routing.overland_flow.parameters - for i in eachindex(routing.overland_flow.storage) - if river[i] + (; bankfull_storage) = routing.river_flow.parameters + for i in eachindex(land_v.storage) + if river_location[i] j = domain.land.network.river_indices[i] if land_v.h[i] > 0.0 land_v.storage[i] = @@ -181,8 +181,7 @@ function set_states!(model::AbstractModel{<:Union{SbmModel, SbmGwfModel}}) land_v.storage[i] = river_v.h[j] * flow_width[j] * flow_length[j] end else - routing.overland_flow.storage[i] = - routing.overland_flow.h[i] * x_length[i] * y_length[i] + land_v.storage[i] = land_v.h[i] * x_length[i] * y_length[i] end end end From 8b9c70c754f97fb17513174fb0a1799722500cbb Mon Sep 17 00:00:00 2001 From: Willem van Verseveld Date: Thu, 22 May 2025 15:07:10 +0200 Subject: [PATCH 2/5] Update and fix local inertial model river and land routing Water depth at waterbody locations (boundary locations) was not fixed at zero during the simulation. Improved the order of computations for the local inertial model river and land routing. --- Wflow/src/routing/surface_local_inertial.jl | 188 ++++++++++++-------- 1 file changed, 115 insertions(+), 73 deletions(-) diff --git a/Wflow/src/routing/surface_local_inertial.jl b/Wflow/src/routing/surface_local_inertial.jl index 6583b0d32..de93c5bbd 100644 --- a/Wflow/src/routing/surface_local_inertial.jl +++ b/Wflow/src/routing/surface_local_inertial.jl @@ -790,8 +790,11 @@ function update!( if t + dt_s > dt dt_s = dt - t end + local_inertial_update_fluxes!(land, domain, dt_s) + update_inflow_waterbody!(land, river, domain) local_inertial_river_update!(river, domain, dt_s, dt, doy, update_h) - local_inertial_update!(land, river, domain, dt_s) + local_inertial_update_water_depth!(land, river, domain, dt_s) + t = t + dt_s end average_flow_vars!(river.variables, actual_external_abstraction_av, dt) @@ -803,27 +806,16 @@ function update!( end """ -Update combined river `LocalInertialRiverFlow`and overland flow `LocalInertialOverlandFlow` -models for a single timestep `dt`. +Update fluxes for overland flow `LocalInertialOverlandFlow` model for a single timestep +`dt`. """ -function local_inertial_update!( +function local_inertial_update_fluxes!( land::LocalInertialOverlandFlow, - river::LocalInertialRiverFlow, domain::Domain, dt::Float64, ) indices = domain.land.network.edge_indices - inds_river = domain.land.network.river_indices - - (; flow_width, flow_length) = domain.river.parameters (; x_length, y_length) = domain.land.parameters - - (; edges_at_node) = domain.river.network - - river_bc = river.boundary_conditions - river_v = river.variables - river_p = river.parameters - land_bc = land.boundary_conditions land_v = land.variables land_p = land.parameters @@ -913,71 +905,121 @@ function local_inertial_update!( end end end + return nothing +end + +""" +Update boundary condition inflow to a waterbody from land `inflow_waterbody` of combined +river `LocalInertialRiverFlow`and overland flow `LocalInertialOverlandFlow` models for a +single timestep. +""" +function update_inflow_waterbody!( + land::LocalInertialOverlandFlow, + river::LocalInertialRiverFlow, + domain::Domain, +) + indices = domain.land.network.edge_indices + inds_river = domain.land.network.river_indices + (; river_location, waterbody_outlet) = domain.land.parameters + + river_bc = river.boundary_conditions + land_bc = land.boundary_conditions + land_v = land.variables + land_p = land.parameters + + @batch per = thread minbatch = 6000 for i in 1:(land_p.n) + yd = indices.yd[i] + xd = indices.xd[i] + if river_location[i] && waterbody_outlet[i] + river_bc.inflow_waterbody[inds_river[i]] = + land_bc.inflow_waterbody[i] + + land_bc.runoff[i] + + (land_v.qx[xd] - land_v.qx[i] + land_v.qy[yd] - land_v.qy[i]) + end + end + return nothing +end + +""" +Update storage and water depth for combined river `LocalInertialRiverFlow`and overland flow +`LocalInertialOverlandFlow` models for a single timestep `dt`. +""" +function local_inertial_update_water_depth!( + land::LocalInertialOverlandFlow, + river::LocalInertialRiverFlow, + domain::Domain, + dt::Float64, +) + indices = domain.land.network.edge_indices + inds_river = domain.land.network.river_indices + (; edges_at_node) = domain.river.network + (; river_location, waterbody_outlet, x_length, y_length) = domain.land.parameters + (; flow_width, flow_length) = domain.river.parameters - # change in storage and water levels based on horizontal fluxes for river and land cells + river_bc = river.boundary_conditions + river_v = river.variables + river_p = river.parameters + land_bc = land.boundary_conditions + land_v = land.variables + land_p = land.parameters + + # change in storage and water depth based on horizontal fluxes for river and land cells @batch per = thread minbatch = 6000 for i in 1:(land_p.n) yd = indices.yd[i] xd = indices.xd[i] - if domain.land.parameters.river_location[i] - if domain.river.parameters.waterbody_outlet[inds_river[i]] - # for reservoir or lake set inflow from land part, these are boundary points - # and update of storage and h is not required - river_bc.inflow_waterbody[inds_river[i]] = - land_bc.inflow_waterbody[i] + - land_bc.runoff[i] + - (land_v.qx[xd] - land_v.qx[i] + land_v.qy[yd] - land_v.qy[i]) + if river_location[i] + # waterbody locations are boundary points (update storage and h not required) + waterbody_outlet[i] && continue + # internal abstraction (water demand) is limited by river storage and negative + # external inflow as part of water allocation computations. + land_v.storage[i] += + ( + sum_at(river_v.q, edges_at_node.src[inds_river[i]]) - + sum_at(river_v.q, edges_at_node.dst[inds_river[i]]) + land_v.qx[xd] - + land_v.qx[i] + land_v.qy[yd] - land_v.qy[i] + land_bc.runoff[i] - + river_bc.abstraction[inds_river[i]] + ) * dt + if land_v.storage[i] < 0.0 + land_v.error[i] = land_v.error[i] + abs(land_v.storage[i]) + land_v.storage[i] = 0.0 # set storage to zero + end + # limit negative external inflow + if river_bc.inflow[inds_river[i]] < 0.0 + available_volume = + if land_v.storage[i] >= river_p.bankfull_storage[inds_river[i]] + river_p.bankfull_depth[inds_river[i]] + else + river_v.storage[inds_river[i]] + end + _inflow = + max(-(0.80 * available_volume / dt), river_bc.inflow[inds_river[i]]) + river_bc.actual_external_abstraction_av[inds_river[i]] += _inflow * dt else - # internal abstraction (water demand) is limited by river storage and negative - # external inflow as part of water allocation computations. - land_v.storage[i] += - ( - sum_at(river_v.q, edges_at_node.src[inds_river[i]]) - - sum_at(river_v.q, edges_at_node.dst[inds_river[i]]) + - land_v.qx[xd] - land_v.qx[i] + land_v.qy[yd] - land_v.qy[i] + - land_bc.runoff[i] - river_bc.abstraction[inds_river[i]] - ) * dt - if land_v.storage[i] < 0.0 - land_v.error[i] = land_v.error[i] + abs(land_v.storage[i]) - land_v.storage[i] = 0.0 # set storage to zero - end - # limit negative external inflow - if river_bc.inflow[inds_river[i]] < 0.0 - available_volume = - if land_v.storage[i] >= river_p.bankfull_storage[inds_river[i]] - river_p.bankfull_depth[inds_river[i]] - else - river_v.storage[inds_river[i]] - end - _inflow = - max(-(0.80 * available_volume / dt), river_bc.inflow[inds_river[i]]) - river_bc.actual_external_abstraction_av[inds_river[i]] += _inflow * dt - else - _inflow = river_bc.inflow[inds_river[i]] - end - land_v.storage[i] += _inflow * dt - if land_v.storage[i] >= river_p.bankfull_storage[inds_river[i]] - river_v.h[inds_river[i]] = - river_p.bankfull_depth[inds_river[i]] + - (land_v.storage[i] - river_p.bankfull_storage[inds_river[i]]) / - (x_length[i] * y_length[i]) - land_v.h[i] = - river_v.h[inds_river[i]] - river_p.bankfull_depth[inds_river[i]] - river_v.storage[inds_river[i]] = - river_v.h[inds_river[i]] * - flow_length[inds_river[i]] * - flow_width[inds_river[i]] - else - river_v.h[inds_river[i]] = - land_v.storage[i] / - (flow_length[inds_river[i]] * flow_width[inds_river[i]]) - land_v.h[i] = 0.0 - river_v.storage[inds_river[i]] = land_v.storage[i] - end - # average variables (here accumulated for model timestep Δt) - river_v.h_av[inds_river[i]] += river_v.h[inds_river[i]] * dt - river_v.storage_av[inds_river[i]] += river_v.storage[inds_river[i]] * dt + _inflow = river_bc.inflow[inds_river[i]] + end + land_v.storage[i] += _inflow * dt + if land_v.storage[i] >= river_p.bankfull_storage[inds_river[i]] + river_v.h[inds_river[i]] = + river_p.bankfull_depth[inds_river[i]] + + (land_v.storage[i] - river_p.bankfull_storage[inds_river[i]]) / + (x_length[i] * y_length[i]) + land_v.h[i] = + river_v.h[inds_river[i]] - river_p.bankfull_depth[inds_river[i]] + river_v.storage[inds_river[i]] = + river_v.h[inds_river[i]] * + flow_length[inds_river[i]] * + flow_width[inds_river[i]] + else + river_v.h[inds_river[i]] = + land_v.storage[i] / + (flow_length[inds_river[i]] * flow_width[inds_river[i]]) + land_v.h[i] = 0.0 + river_v.storage[inds_river[i]] = land_v.storage[i] end + # average variables (here accumulated for model timestep Δt) + river_v.h_av[inds_river[i]] += river_v.h[inds_river[i]] * dt + river_v.storage_av[inds_river[i]] += river_v.storage[inds_river[i]] * dt else land_v.storage[i] += ( From b5dcd48092f421b6c48afdfd11b4868352d48b5c Mon Sep 17 00:00:00 2001 From: Willem van Verseveld Date: Thu, 22 May 2025 15:37:38 +0200 Subject: [PATCH 3/5] Update changelog --- docs/changelog.qmd | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/changelog.qmd b/docs/changelog.qmd index 2f924cf92..876cd17dc 100644 --- a/docs/changelog.qmd +++ b/docs/changelog.qmd @@ -31,6 +31,8 @@ project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). flow) for river cells in the diagonal flow direction. For both cases the width should be based on land area and flow length (area conservative). - Limit external negative inflow for local inertial routing (maximum of 80% of river volume). +- Water depth for waterbody locations (boundary condition) of the local inertial model for + river and land routing was computed and should remain zero during the simulation. ### Changed - Removed vertical concepts `HBV` and `FLEXTopo`. From e1b450f7424700bbfaa746e9c22a81364071a916 Mon Sep 17 00:00:00 2001 From: Willem van Verseveld Date: Mon, 26 May 2025 13:05:26 +0200 Subject: [PATCH 4/5] Update changelog --- docs/changelog.qmd | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/docs/changelog.qmd b/docs/changelog.qmd index 876cd17dc..91e4072d0 100644 --- a/docs/changelog.qmd +++ b/docs/changelog.qmd @@ -31,8 +31,6 @@ project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). flow) for river cells in the diagonal flow direction. For both cases the width should be based on land area and flow length (area conservative). - Limit external negative inflow for local inertial routing (maximum of 80% of river volume). -- Water depth for waterbody locations (boundary condition) of the local inertial model for - river and land routing was computed and should remain zero during the simulation. ### Changed - Removed vertical concepts `HBV` and `FLEXTopo`. @@ -103,6 +101,14 @@ project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). For the transpiration from the saturated store the remaining potential transpiration is used in combination with the wet root fraction (root length density fraction in the saturated soil layers is not used). +- Changed order of computations for 1D-2D local inertial model that improves the accuracy of + reservoir or lake state after initialization with a warm state. The inflow from land to a + reservoir or lake was updated during the overland flow routing (final computation during + each internal time step). This results in small inflow deviations when starting with a + warm state as the inflow from land to a reservoir or lake has not been updated yet, while + a reservoir or lake is updated as part of the river routing. Order of computations has + been changed to: 1) overland flow routing, 2) reservoir or lake inflow and 3) river + routing. ### Added - Support direct output of snow and glacier melt, and add computation of snow water From 4f3208e5717827c8340182f23f43575014f1ad36 Mon Sep 17 00:00:00 2001 From: Willem van Verseveld Date: Mon, 26 May 2025 14:02:17 +0200 Subject: [PATCH 5/5] Consistent use of `check_timestepsize` function --- Wflow/src/routing/surface_local_inertial.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/Wflow/src/routing/surface_local_inertial.jl b/Wflow/src/routing/surface_local_inertial.jl index de93c5bbd..535c3a3dd 100644 --- a/Wflow/src/routing/surface_local_inertial.jl +++ b/Wflow/src/routing/surface_local_inertial.jl @@ -500,9 +500,7 @@ function update!( t = 0.0 while t < dt dt_s = stable_timestep(model, flow_length) - if t + dt_s > dt - dt_s = dt - t - end + dt_s = check_timestepsize(dt_s, t, dt) local_inertial_river_update!(model, domain, dt_s, dt, doy, update_h) t = t + dt_s end @@ -787,9 +785,8 @@ function update!( dt_river = stable_timestep(river, flow_length) dt_land = stable_timestep(land, parameters) dt_s = min(dt_river, dt_land) - if t + dt_s > dt - dt_s = dt - t - end + dt_s = check_timestepsize(dt_s, t, dt) + local_inertial_update_fluxes!(land, domain, dt_s) update_inflow_waterbody!(land, river, domain) local_inertial_river_update!(river, domain, dt_s, dt, doy, update_h)