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
+    reinfiltration_surfacewater__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,35 @@ 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 correct_overland_flow_level!(
+    model::SbmSoilModel,
+    overland_flow::KinWaveOverlandFlowModel,
+    domain::Domain,
+    config::Config,
+)
+    v = model.variables
+
+    if config.model.reinfiltration_surfacewater__flag
+        (; surface_flow_width) = domain.land.parameters
+        n = length(surface_flow_width)
+        threaded_foreach(1:n; basesize = 1000) do i
+            q, h = correct_overland_flow_level(
+                overland_flow.variables.h[i],
+                v.infilt_surfacewater[i],
+                domain.land.parameters.river_fraction[i],
+                surface_flow_width[i],
+                overland_flow.parameters.alpha[i],
+                overland_flow.parameters.beta,
+            )
+            if !isnothing(q)
+                overland_flow.variables.flow.q[i] = q
+                overland_flow.variables.h[i] = h
+            end
+        end
+    end
+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"
@@ -1020,6 +1021,11 @@ function local_inertial_update_water_depth!(
                     land_v.qx[i] + land_v.qy[yd] - land_v.qy[i] + land_bc.runoff[i] -
                     river_bc.abstraction[inds_river[i]]
                 ) * dt
+            # Apply surface water infiltration correction for river cells
+            if land_bc.infiltration_volume[i] > 0.0
+                land_v.storage[i] =
+                    max(0.0, land_v.storage[i] - land_bc.infiltration_volume[i])
+            end
             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
@@ -1070,6 +1076,11 @@ function local_inertial_update_water_depth!(
                 land_v.error[i] = land_v.error[i] + abs(land_v.storage[i])
                 land_v.storage[i] = 0.0 # set storage to zero
             end
+            # Apply surface water infiltration correction if available
+            if land_bc.infiltration_volume[i] > 0.0
+                land_v.storage[i] =
+                    max(0.0, land_v.storage[i] - land_bc.infiltration_volume[i])
+            end
             land_v.h[i] = land_v.storage[i] / (x_length[i] * y_length[i])
         end
     end
@@ -1353,3 +1364,25 @@ function FloodPlainModel(
     floodplain = FloodPlainModel(; parameters, variables)
     return floodplain
 end
+
+"""
+Update overland flow water level and storage for LocalInertialOverlandFlow model based on
+surface water infiltration. For local inertial flow, the infiltration is applied by updating
+the boundary conditions which are then used in the water depth update function.
+"""
+function correct_overland_flow_level!(
+    model::SbmSoilModel,
+    overland_flow::LocalInertialOverlandFlowModel,
+    domain::Domain,
+    config::Config,
+)
+    (; infilt_surfacewater) = model.variables
+    (; area) = domain.land.parameters
+
+    if config.model.reinfiltration_surfacewater__flag
+        # Update the boundary condition for surface water infiltration
+        # This will be used in local_inertial_update_water_depth!
+        overland_flow.boundary_conditions.infiltration_volume .=
+            infilt_surfacewater .* area .* 0.001
+    end
+end

Wflow/src/sbm.jl

@@ -125,7 +125,15 @@ function update_land_hydrology_model!(
         (; interception, runoff, demand, allocation),
     )
 
-    update_soil_water_flow!(soil, atmospheric_forcing, (; snow, runoff, demand), config, dt)
+    update_soil_water_flow!(
+        soil,
+        domain,
+        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,12 @@ 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)
+    # 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)
     # 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⁻¹]
@@ -687,12 +693,13 @@ function update_bc_soil_model!(
     evaporation!(demand.paddy, potential_soilevaporation)
     potential_soilevaporation .= potential_soilevaporation .- get_evaporation(demand.paddy)
 
-    water_flux_surface .= max.(
-        runoff.boundary_conditions.water_flux_surface .+
-        get_irrigation_allocated(allocation) .- runoff.variables.runoff_river .-
-        runoff.variables.runoff_land .+ get_water_depth(demand.paddy),
-        0.0,
-    )
+    water_flux_surface .=
+        max.(
+            runoff.boundary_conditions.water_flux_surface .+
+            get_irrigation_allocated(allocation) .- runoff.variables.runoff_river .-
+            runoff.variables.runoff_land .+ get_water_depth(demand.paddy),
+            0.0,
+        )
     return nothing
 end
 
@@ -741,6 +748,51 @@ 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) = model.boundary_conditions
+    (; waterdepth_land) = runoff.boundary_conditions
+    (; river_fraction) = domain.land.parameters
+
+    n = length(v.potential_infiltration)
+    threaded_foreach(1:n; basesize = 1000) do i
+        v.potential_infiltration_surfacewater[i] = 0.0
+        if do_surface_water_infiltration
+            v.potential_infiltration_surfacewater[i] =
+                waterdepth_land[i] * (1.0 - river_fraction[i]) * 0.95
+            water_flux_surface[i] += v.potential_infiltration_surfacewater[i]
+        end
+        v.potential_infiltration[i] = water_flux_surface[i]
+    end
+
+    return nothing
+end
+
+function correct_infiltration!(model::SbmSoilModel)
+    v = model.variables
+    (; water_flux_surface) = model.boundary_conditions
+
+    n = length(v.actinfilt)
+    threaded_foreach(1:n; basesize = 1000) do i
+        v.infilt_surfacewater[i],
+        v.actinfilt[i],
+        v.infiltexcess[i],
+        v.excesswater[i],
+        water_flux_surface[i] = correct_infiltration(
+            v.potential_infiltration[i],
+            v.potential_infiltration_surfacewater[i],
+            water_flux_surface[i],
+            v.actinfilt[i],
+            v.infiltexcess[i],
+        )
+    end
+end
+
 """
     infiltration!(soil_model::SbmSoilMsoil
 
@@ -1118,6 +1170,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,
@@ -1138,6 +1191,14 @@ function update_soil_water_flow!(
         modelsnow = config.model.snow__flag,
         soil_infiltration_reduction = config.model.soil_infiltration_reduction__flag,
     )
+    # update available for infiltration in case surface water infiltration is enabled
+    update_available_for_infiltration!(
+        soil_model,
+        domain,
+        runoff,
+        config.model.reinfiltration_surfacewater__flag,
+    )
+
     infiltration!(soil_model)
     # unsaturated zone flow
     unsaturated_zone_flow!(soil_model)
@@ -1146,8 +1207,13 @@ 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, also to ensure correct soil and path
+    # infiltration, and excesswater
+    correct_infiltration!(soil_model)
+
     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 +1304,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 +1391,10 @@ 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
+
+    # correct overland flow water levels in case of reinfiltration
+    correct_overland_flow_level!(soil_model, overland_flow, domain, config)
+
     return nothing
 end
 

Wflow/src/soil/soil_process.jl

@@ -267,3 +267,70 @@ function actual_infiltration_soil_path(
 
     return actinfiltsoil, actinfiltpath
 end
+
+"""
+Correct infiltration fluxes by separating the surface water contribution from the total
+infiltration. The correction factor is based on the ratio of
+`potential_infiltration_surfacewater` to `potential_infiltration`, and is applied to
+`actinfilt` and `infiltexcess` to remove the surface water component. The surface water flux
+`water_flux_surface` is adjusted accordingly, and the remaining excess water is computed.
+
+Returns `infilt_surfacewater`, corrected `actinfilt`, corrected `infiltexcess`, `excesswater`,
+and corrected `water_flux_surface`.
+"""
+function correct_infiltration(
+    potential_infiltration,
+    potential_infiltration_surfacewater,
+    water_flux_surface,
+    actinfilt,
+    infiltexcess,
+)
+    # Determine ratio of water that has infiltrated
+    infilt_ratio = potential_infiltration == 0.0 ? 0.0 : actinfilt / potential_infiltration
+    # Use this ratio to determine the contribution from overland flow
+    infilt_surfacewater = max(0.0, potential_infiltration_surfacewater * infilt_ratio)
+    # Determine the correction factor to apply to the relevant fluxes
+    correction_surfacewater =
+        potential_infiltration == 0.0 ? 1.0 :
+        1.0 - (potential_infiltration_surfacewater / potential_infiltration)
+
+    # Correct fluxes
+    actinfilt *= correction_surfacewater
+    infiltexcess *= correction_surfacewater
+
+    # subtract contribution from overland flow to ensure correct fluxes
+    water_flux_surface -= potential_infiltration_surfacewater
+    excesswater = water_flux_surface - actinfilt - infiltexcess
+
+    return infilt_surfacewater, actinfilt, infiltexcess, excesswater, water_flux_surface
+end
+
+function correct_overland_flow_level(
+    overlandflow_depth,
+    infilt_surfacewater,
+    river_fraction,
+    surface_flow_width,
+    alpha,
+    beta,
+)
+    if infilt_surfacewater > 0.0
+        # Get original h_land in mm
+        original_h_land = overlandflow_depth * 1000.0
+
+        # Correct values for river fraction to ensure correct water accounting
+        infiltrated_surfacewater = (infilt_surfacewater / (1.0 - river_fraction))
+        # Calculate new h_land in m
+        h = (original_h_land - infiltrated_surfacewater) / 1000.0
+
+        q = ifelse(
+            surface_flow_width > 0.0 && h > 0.0 && alpha > 0.0 && beta != 0.0,
+            # Compute cross-sectional area from h
+            pow.(max.((h * surface_flow_width) / alpha, 1e-10), 1.0 / beta),
+            0.0,  # Set q to 0.0 if conditions are not met
+        )
+    else
+        q = nothing
+        h = overlandflow_depth
+    end
+    return q, h
+end