Decide on the hydrological concepts of the agricultural module

[deltamarnix]

Which wflow concepts need to go into the shared library? What does Ribasim need? What does an agricultural library need to function?

• Precipitation
• Interception
• Soil column boundary conditions

[vers-w]

The concepts/structs that are at least required are the following (from LandHydrologySBM):

• atmospheric_forcing
• vegetation_parameters
• interception
• soil

The variable water_table_depth is a boundary condition for unsaturated flow (set by subsurface flow). The Ribasim agricultural module most probably requires another boundary condition (e.g. free drainage). Some processes involving the saturated zone as capillary rise and transpiration from the saturated store (wet roots) are probably also not needed.

[HendrikKok]

The shared irrigation module between Wflow and Ribasim could have the following setup. We define a 1D SBM soil column per irrigation water user in Ribasim. Each water user is represented by a single-column LandHydrologySBM model with one crop type. This model has the following setup:

Active:

1. GashInterceptionModel — canopy interception, throughfall, stemflow
2. SbmSoilModel — infiltration, unsaturated zone drainage (Brooks-Corey), soil evaporation, transpiration (Feddes stress), capillary rise, leakage to deep groundwater, water table dynamics
3. OpenWaterRunoff — passes throughfall to the soil surface boundary condition
4. DemandModel (NonPaddy) — computes soil-moisture-deficit-based irrigation demand internally
5. ExternalIrrigationAllocation — injects externally allocated irrigation flux into the soil surface

Inactive:

1. NoSnowModel — no snow accumulation or melt; also disables soil temperature tracking
2. NoGlacierModel — no glacier melt
3. NoNonIrrigationDemandModel — no sectoral (domestic/industry/livestock) demand
4. NoSubsurfaceFlow (stub) — no lateral groundwater exchange between cells

Disabled flags:

• snow__flag = false — soil temperature update skipped
• soil_infiltration_reduction__flag = false — infiltration is never reduced for frozen soil
• leaf_area_index = nothing — static canopy, no seasonal LAI variation

The column can have 1 to n layers, but always represents a single crop type.

Time-dependent inputs:

1. Precipitation
2. Potential evapotranspiration (PET)
3. Temperature (can be constant; affects infiltration reduction for frozen soil if enabled)
4. Leaf area index / LAI (optional, for seasonal canopy dynamics)

Per-layer or column-level parameters:

1. Vegetation/crop — rooting depth, root fraction per layer, crop coefficient Kc
2. Soil hydraulics — Ksat at surface, decay scale, θs, θr, θfc, Brooks-Corey exponent, infiltration capacity, maximum leakage
3. Feddes stress function — h1, h2, h3_high, h3_low, h4

This column will be the core of the irrigation module and will be wrapped together with some additional methods.

Methods:

1. initialise(params) — constructs the SBM model from crop and soil parameters
2. get_demand() — calls update_water_demand_model! on the current soil state θ(t-1), returns demand_gross [m/s]
3. set_allocated(flux) — sets allocation.irrigation_allocation = flux [m/s]; must be called before update_until
4. update_until(t, dt) — advances the soil from t to t+dt:
   • Reset allocation.irrigation_allocation = 0.0
   • Set atmospheric_forcing for t
   • update_interception_model!
   • get_water_flux_surface!
   • update_bc_soil_model! ← allocated irrigation enters the soil here
   • update_soil_water_flow!
   • update_water_table!
   • update_soil_water_storage!
5. update_water_table! — new method replacing kinematic_wave_ssf / update_subsurface_flow_model!. Uses water_table_change with specific_yield = lower_bound_drainable_porosity(θs, θfc) (not θs − θr).
6. get_runoff() — returns net_runoff [m/s]

Additional thoughts: irrigation module keeps track of its own time. Time dependent input could be provided at initialization, or could be set externally per exchange step.

External coupling sequence:

irrigation.initialise(params)

# Per exchange timestep Δt
demand    = irrigation.get_demand()       # based on θ(t-1)
irrigation.set_allocated(allocated)       # must be called before update_until
irrigation.update_until(t+dt)           # advances soil from t to t + Δt

# Optional: return fluxes to Ribasim
irrigation.get_runoff()

t += dt

Example output (single-column test, 30-day run):

Top to bottom: daily water budget (inputs/outputs balanced); water table depth;
irrigation demand vs allocation (50% supply fraction); cumulative fluxes.

Note on timing: get_demand() and set_allocated() must be called before update_until. A natural one-step lag applies: demand is based on θ(t-1) and supply arrives within the same timestep. A full one-step lag can be implemented by deferring set_allocated to the following timestep.

Potential extensions:

• Route get_runoff() to basins.inflow in Ribasim

The module uses Wflow internals and adds a small number of additional methods. Wflow itself has no dependency on this module.

[gijsber]

Nice write up @HendrikKok
On the Potential extension:
Shouldn't get_runoff() go to basins.surface_runoff in Ribasim instead ?

Does the call stack account for the reduced crop survival when the there is a soil water defiicit ? (if so, which call is doing it ?)

[deltamarnix]

For a proof of concept I can live with the irrigation module referencing to wflow, but from an architecture viewpoint it makes more sense to make the irrigation module standalone, so that both Ribasim and Wflow can point to the same module.

One of the investigations we should do is to see how we can separate out the irrigation module as a standalone entity. Otherwise Ribasim completely includes Wflow as a reference just for a single module. I can imagine that the setup would also become more cumbersome for a Ribasim model.

[vers-w]

I agree that as a first step for the prototype it makes sense to make no changes on the Wflow side: Deltares/Wflow.jl#935 (comment). However, when new functionality is required (e.g. using dynamic crop factors, or different crop types per grid cell/irrigation node) it would be good to add this to Wflow so this is shared.

[HendrikKok]

@gijsber: Yes, that is indeed the correct address.

@deltamarnix:
So far, I have mainly focused on identifying which concepts from Wflow are needed for this module. The current approach is to use the irrigation module as a wrapper around the SBM model concepts. Within Wflow, the SBM model is used in a more complex setup, so I do not see a need for Wflow itself to use the SBM concepts through the irrigation wrapper. As @vers-w mentioned, additions to the SBM core should indeed be implemented in the Wflow codebase.

An alternative would be to extract the SBM model into a separate shared component. Both Wflow and the irrigation module could then depend on the same implementation, with updates to the SBM core being made in that shared codebase. This would provide a clearer separation between the core model and its applications I guess?

[HendrikKok]

So the first draft of the prototype is located here: #3148. It includes some stand-alone examples on the api. Together with @SouthEndMusic we will try to run this online with Ribasim.

[HendrikKok]

I had a look at the hydrological concepts in the SBM-model, regarding the scope of this module. So from top down we have the following potential issues:

Bare Soil Evaporation (code change)
The current approach scales potential evaporation by the relative wetness of the top layer, making results discretisation-dependent. Suggest adopting Boesten & Stroosnijder (1986) as used in SWAP/MetaSWAP.

ET from Root Zone (config/parameter)
Distributing transpiration over unsaturated layers by root fraction is sound. Direct extraction from the saturated zone is no. The capillary fringe is under tension and should be treated as unsaturated storage. The correct pathway is:

capillary rise → unsaturated storage → root uptake

Set soil_wet_root__sigmoid_function_shape_parameter to a large value to suppress saturated-zone transpiration.

Percolation (config)
The exponential ksat decay with depth lacks a physical basis. Always use the layered ksat profile, even for a single layer, to bypass this.

Inter-layer Flux (no change needed)
k(θ) with unit head gradient is acceptable for gravity-driven drainage.

Lower Boundary (code change)
Current implementation uses ksat with unit gradient throughout. When the water table is below the soil column, k(θ) should be used instead of ksat. When the water table is inside the column, ksat(=k(θ)) is appropriate but the unit gradient is not. Suggest introducing a Cauchy (General Head) boundary driven by an external aquifer head.

Capillary Rise (code change)
Current implementation uses ksat at the water table, which heavily overestimates upward flux. Suggest using k(θ) with the actual pressure head gradient between the saturated zone and the bottom unsaturated layer.

[deltamarnix]

@HendrikKok Can you create a separate issue on these final points and close this issue?

[HendrikKok]

Created Deltares/Wflow.jl#995 to further look into the code change issues