Surfacewater infiltration

Wflow/src/config_structure.jl

@@ -70,6 +70,7 @@ end
     snow_gravitational_transport__flag::Bool = false
     glacier__flag::Bool = false
     soil_infiltration_reduction__flag::Bool = false
+    land_surface_water_reinfiltration__flag::Bool = false
     soil_layer__thickness::Vector{Int} = [100, 300, 800]
     saturated_hydraulic_conductivity_profile::VerticalConductivityProfile.T =
         VerticalConductivityProfile.exponential

Wflow/src/routing/surface/surface_kinwave.jl

@@ -649,3 +649,63 @@ get_inflow_reservoir(
 # Exclude subsurface flow from `GroundwaterFlowModel`.
 get_inflow_reservoir(::AbstractRiverFlowModel, ::GroundwaterFlowModel, inds::Vector{Int}) =
     zeros(length(inds))
+
+"""
+Update overland flow water level and discharge for KinWaveOverlandFlow model based on
+surface water infiltration.
+"""
+function update_overland_flow_and_depth!(
+    overland_flow_model::KinWaveOverlandFlowModel,
+    soil_model::SbmSoilModel,
+    domain::Domain,
+)
+    (; infilt_surfacewater) = soil_model.variables
+    n = length(infilt_surfacewater)
+    threaded_foreach(1:n; basesize = 1000) do i
+        update_overland_flow_and_depth!(
+            overland_flow_model,
+            infilt_surfacewater[i],
+            domain.land.parameters,
+            i,
+        )
+    end
+end
+
+"""
+Update overland flow water level and discharge in-place for a single cell based on
+surface water infiltration.
+"""
+function update_overland_flow_and_depth!(
+    overland_flow_model::KinWaveOverlandFlowModel,
+    infilt_surfacewater,
+    land_parameters,
+    i,
+)
+    if infilt_surfacewater > 0.0
+        # Get original h_land in mm
+        original_h_land = overland_flow_model.variables.h[i] * 1000.0
+
+        # Correct values for river fraction to ensure correct water accounting
+        infiltrated_surfacewater =
+            (infilt_surfacewater / (1.0 - land_parameters.river_fraction[i]))
+        # Calculate new h_land in m
+        h = (original_h_land - infiltrated_surfacewater) / 1000.0
+
+        q = ifelse(
+            land_parameters.surface_flow_width[i] > 0.0 && h > 0.0,
+            # Compute cross-sectional area from h
+            pow(
+                (h * land_parameters.surface_flow_width[i]) /
+                overland_flow_model.parameters.alpha[i],
+                1.0 / overland_flow_model.parameters.beta,
+            ),
+            0.0,
+        )
+        # set q to 0.0 if it is below the minimum flow threshold
+        q = ifelse(q < KIN_WAVE_MIN_FLOW, 0.0, q)
+
+        overland_flow_model.variables.flow.q[i] = q
+        overland_flow_model.variables.h[i] = h
+    end
+    return nothing
+end

Wflow/src/routing/surface/surface_local_inertial.jl

@@ -655,6 +655,7 @@ end
 @with_kw struct LocalInertialOverlandFlowBC
     n::Int
     runoff::Vector{Float64} = zeros(n) # runoff from hydrological model [m³ s⁻¹]
+    infiltration_volume::Vector{Float64} = zeros(n) # amount of infiltration from surface water [m³]
 end
 
 "Local inertial overland flow model using the local inertial method"
@@ -752,6 +753,10 @@ function update_bc_overland_flow_model!(
         net_runoff / 1000.0 * area / dt + net_runoff_river * area * 0.001 / dt
     overland_flow_model.boundary_conditions.runoff[river_indices] .+=
         get_flux_to_river(subsurface_flow, river_indices)
+
+    # infiltration volume [m³] from surface water, applied in `local_inertial_update_water_depth!`
+    @. overland_flow_model.boundary_conditions.infiltration_volume =
+        soil.variables.infilt_surfacewater * area * 0.001
     return nothing
 end
 
@@ -1019,7 +1024,7 @@ function local_inertial_update_water_depth!(
                     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
+                ) * dt - land_bc.infiltration_volume[i]
             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
@@ -1065,7 +1070,7 @@ function local_inertial_update_water_depth!(
                 (
                     land_v.qx[xd] - land_v.qx[i] + land_v.qy[yd] - land_v.qy[i] +
                     land_bc.runoff[i]
-                ) * dt
+                ) * dt - land_bc.infiltration_volume[i]
             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
@@ -1353,3 +1358,14 @@ function FloodPlainModel(
     floodplain = FloodPlainModel(; parameters, variables)
     return floodplain
 end
+
+"""
+No-op for the `LocalInertialOverlandFlowModel`. Surface water infiltration is applied via
+the `infiltration_volume` boundary condition set in `update_bc_overland_flow_model!` and
+subtracted from land storage in `local_inertial_update_water_depth!`.
+"""
+update_overland_flow_and_depth!(
+    ::LocalInertialOverlandFlowModel,
+    ::SbmSoilModel,
+    ::Domain,
+) = nothing

Wflow/src/routing/surface/surface_routing.jl

@@ -10,6 +10,11 @@ function surface_routing!(model)
     (; overland_flow, river_flow, subsurface_flow) = routing
     (; reservoir) = river_flow.boundary_conditions
 
+    # correct overland flow water levels in case of reinfiltration
+    if config.model.land_surface_water_reinfiltration__flag
+        update_overland_flow_and_depth!(overland_flow, soil, domain)
+    end
+
     dt = tosecond(clock.dt)
     # update lateral inflow for kinematic wave overland flow
     update_lateral_inflow!(

Wflow/src/sbm.jl

@@ -123,9 +123,19 @@ function update_land_hydrology_model!(
         soil,
         atmospheric_forcing,
         (; interception, runoff, demand, allocation),
+        domain,
+        config,
+    )
+
+    update_soil_water_flow!(
+        soil,
+        domain,
+        atmospheric_forcing,
+        (; snow, runoff, demand),
+        config,
+        dt,
     )
 
-    update_soil_water_flow!(soil, atmospheric_forcing, (; snow, runoff, demand), config, dt)
     @. soil.variables.actevap += interception.variables.interception_rate
     return nothing
 end

Wflow/src/sbm_gwf_model.jl

@@ -64,7 +64,8 @@ end
 function update_model!(model::AbstractModel{<:SbmGwfModel})
     (; routing, land, domain, clock, config) = model
     (; soil, runoff, demand) = land
-    (; boundary_conditions) = routing.subsurface_flow
+    (; subsurface_flow, overland_flow) = routing
+    (; boundary_conditions) = subsurface_flow
 
     dt = tosecond(clock.dt)
 
@@ -85,7 +86,7 @@ function update_model!(model::AbstractModel{<:SbmGwfModel})
     end
     # update groundwater domain
     update_subsurface_flow_model!(
-        routing.subsurface_flow,
+        subsurface_flow,
         soil,
         domain,
         dt_gwf,
@@ -94,7 +95,9 @@ function update_model!(model::AbstractModel{<:SbmGwfModel})
     # update SBM soil model (runoff, ustorelayerdepth and satwaterdepth)
     update_soil_water_storage!(
         soil,
-        (; runoff, demand, subsurface_flow = routing.subsurface_flow),
+        (; runoff, demand, subsurface_flow, overland_flow),
+        domain,
+        config,
     )
 
     surface_routing!(model)

Wflow/src/sbm_model.jl

@@ -61,6 +61,7 @@ function update_model!(model::AbstractModel{<:SbmModel})
     (; routing, land, domain, clock, config) = model
     (; soil, runoff, demand) = land
     (; kv_profile) = land.soil.parameters
+    (; subsurface_flow, overland_flow) = routing
     (; boundary_conditions) = routing.subsurface_flow
     dt = tosecond(clock.dt)
 
@@ -85,8 +86,14 @@ function update_model!(model::AbstractModel{<:SbmModel})
         domain,
         clock.dt / BASETIMESTEP,
     )
+
     # update SBM soil model (runoff, ustorelayerdepth and satwaterdepth)
-    update_soil_water_storage!(soil, (; runoff, demand, routing.subsurface_flow))
+    update_soil_water_storage!(
+        soil,
+        (; runoff, demand, subsurface_flow, overland_flow),
+        domain,
+        config,
+    )
 
     surface_routing!(model)
 

Wflow/src/soil/soil.jl

@@ -43,6 +43,8 @@ abstract type AbstractSoilModel end
     infiltsoilpath::Vector{Float64} = fill(MISSING_VALUE, n)
     # Infiltration excess water [mm Δt⁻¹]
     infiltexcess::Vector{Float64} = fill(MISSING_VALUE, n)
+    # Infiltration from surface water [mm Δt⁻¹]
+    infilt_surfacewater::Vector{Float64} = fill(0.0, n)
     # Water that cannot infiltrate due to saturated soil (saturation excess) [mm Δt⁻¹]
     excesswater::Vector{Float64} = fill(MISSING_VALUE, n)
     # Water exfiltrating during saturation excess conditions [mm Δt⁻¹]
@@ -198,6 +200,10 @@ end
     potential_transpiration::Vector{Float64} = fill(MISSING_VALUE, n)
     # Potential soil evaporation rate [mm Δt⁻¹]
     potential_soilevaporation::Vector{Float64} = fill(MISSING_VALUE, n)
+    # Potential infiltration originating from surface water [mm Δt⁻¹]
+    potential_infiltration_surfacewater::Vector{Float64} = fill(0.0, n)
+    # Total water available for infiltration [mm Δt⁻¹]
+    potential_infiltration::Vector{Float64} = fill(0.0, n)
 end
 
 "Exponential depth profile of vertical hydraulic conductivity at the soil surface"
@@ -675,6 +681,8 @@ function update_bc_soil_model!(
     soil_model::SbmSoilModel,
     atmospheric_forcing::AtmosphericForcing,
     external_models::NamedTuple,
+    domain::Domain,
+    config::Config,
 )
     (; interception, runoff, demand, allocation) = external_models
     (; potential_transpiration, water_flux_surface, potential_soilevaporation) =
@@ -693,6 +701,13 @@ function update_bc_soil_model!(
         runoff.variables.runoff_land .+ get_water_depth(demand.paddy),
         0.0,
     )
+    # update available for infiltration in case surface water infiltration is enabled
+    update_available_for_infiltration!(
+        soil_model,
+        domain,
+        runoff,
+        config.model.land_surface_water_reinfiltration__flag,
+    )
     return nothing
 end
 
@@ -741,6 +756,56 @@ function infiltration_reduction_factor!(
     return nothing
 end
 
+function update_available_for_infiltration!(
+    model::SbmSoilModel,
+    domain::Domain,
+    runoff::AbstractRunoffModel,
+    do_surface_water_infiltration::Bool,
+)
+    v = model.variables
+    (; water_flux_surface, potential_infiltration_surfacewater, potential_infiltration) =
+        model.boundary_conditions
+    (; waterdepth_land) = runoff.boundary_conditions
+    (; river_fraction) = domain.land.parameters
+
+    # Fraction of overland flow that is made available for reinfiltration [-]
+    max_reinfiltration_fraction = 0.95
+
+    n = length(potential_infiltration)
+    threaded_foreach(1:n; basesize = 1000) do i
+        potential_infiltration_surfacewater[i] = 0.0
+        if do_surface_water_infiltration
+            potential_infiltration_surfacewater[i] =
+                waterdepth_land[i] * (1.0 - river_fraction[i]) * max_reinfiltration_fraction
+            water_flux_surface[i] += potential_infiltration_surfacewater[i]
+        end
+        potential_infiltration[i] = water_flux_surface[i]
+    end
+
+    return nothing
+end
+
+function update_infiltration_fluxes!(soil_model::SbmSoilModel)
+    (; infilt_surfacewater, actinfilt, infiltexcess, excesswater) = soil_model.variables
+    (; water_flux_surface, potential_infiltration, potential_infiltration_surfacewater) =
+        soil_model.boundary_conditions
+
+    n = length(actinfilt)
+    threaded_foreach(1:n; basesize = 1000) do i
+        infilt_surfacewater[i],
+        actinfilt[i],
+        infiltexcess[i],
+        excesswater[i],
+        water_flux_surface[i] = update_infiltration_fluxes(
+            potential_infiltration[i],
+            potential_infiltration_surfacewater[i],
+            water_flux_surface[i],
+            actinfilt[i],
+            infiltexcess[i],
+        )
+    end
+end
+
 """
     infiltration!(soil_model::SbmSoilMsoil
 
@@ -1118,6 +1183,7 @@ transpiration, capillary flux and leakage) for a single timestep.
 """
 function update_soil_water_flow!(
     soil_model::SbmSoilModel,
+    domain::Domain,
     atmospheric_forcing::AtmosphericForcing,
     external_models::NamedTuple,
     config::Config,
@@ -1146,8 +1212,17 @@ function update_soil_water_flow!(
     transpiration!(soil_model, dt)
     # actual infiltration and excess water
     actual_infiltration!(soil_model)
-    @. v.excesswater = water_flux_surface - v.actinfilt - v.infiltexcess
+
+    # Correct fluxes in case of reinfiltration, and only compute excess water if reinfiltration
+    # is not enabled
+    if config.model.land_surface_water_reinfiltration__flag
+        update_infiltration_fluxes!(soil_model)
+    else
+        @. v.excesswater = water_flux_surface - v.actinfilt - v.infiltexcess
+    end
+
     actual_infiltration_soil_path!(soil_model)
+
     @. v.excesswatersoil =
         max(water_flux_surface * (1.0 - p.pathfrac) - v.actinfiltsoil, 0.0)
     @. v.excesswaterpath = max(water_flux_surface * p.pathfrac - v.actinfiltpath, 0.0)
@@ -1238,8 +1313,14 @@ store `satwaterdepth` and the water exfiltrating during saturation excess condit
 `exfiltsatwater` are updated. Additionally, volumetric water content per soil layer and for
 the root zone are updated.
 """
-function update_soil_water_storage!(soil_model::SbmSoilModel, external_models::NamedTuple)
-    (; runoff, demand, subsurface_flow) = external_models
+function update_soil_water_storage!(
+    soil_model::SbmSoilModel,
+    external_models::NamedTuple,
+    domain::Domain,
+    config::Config,
+)
+    (; runoff, demand, subsurface_flow, overland_flow) = external_models
+
     (; runoff_land, ae_openw_l) = runoff.variables
     p = soil_model.parameters
     v = soil_model.variables
@@ -1319,6 +1400,7 @@ function update_soil_water_storage!(soil_model::SbmSoilModel, external_models::N
     # and the h_max parameter of a paddy field)
     update_runoff!(demand.paddy, v.runoff)
     @. v.net_runoff = v.runoff - ae_openw_l
+
     return nothing
 end