Skip to content

Irrigation system

An irrigation system, a reservoir and an aquifer that keep the root zone of the crops supplied.

  • Name: irrigation
  • Field: engineering/processes
  • Formalism: ordinary differential equations
  • References: Bala et al. 2017

Factory

@register_model(
    field='engineering/processes',
    tags=('water resources', 'agriculture', 'stocks and flows', 'system dynamics', 'ode'),
    references=('Bala et al. 2017',),
    defaults={'until': 730.0, 'dt': 1.0},
    observables=lambda system: [
        ObservableSpec('demand', IrrigationSchemeElement._demand, element='scheme'),
        ObservableSpec('release', IrrigationSchemeElement._reservoir_release, element='scheme'),
        ObservableSpec('pumping', IrrigationSchemeElement._pumping, element='scheme'),
    ],
)
def irrigation(
    *,
    river_inflow: float = 0.5,
    seasonality: float = 0.9,
    reservoir_capacity: float = 80.0,
    potential_et: float = 0.4,
    canal_loss: float = 0.2,
    pump_capacity: float = 0.2,
) -> System:
    """Return an irrigation system, a reservoir and an aquifer that keep the root zone of the crops supplied.

    A river fills the reservoir with a yearly cycle, and the reservoir loses water to evaporation and spills above its
    capacity. The irrigation demand is the potential crop evapotranspiration plus the refill of the root zone towards
    a target fraction of its field capacity. The reservoir releases the demand grossed up for the canal losses, as far
    as its stock allows, and the losses seep into the aquifer. When the release falls short, the farmers pump the
    unmet demand from the aquifer up to the pump capacity. The crops draw water from the root zone at the potential
    rate scaled by a water-stress table, and the root zone percolates to the aquifer, which drains as baseflow. Water
    is conserved, so in the dry season the reservoir falls, the pumping rises, and the aquifer is drawn down. The
    structure follows the irrigation example of Bala, Arshad and Noh (2017, Section 4.11), with illustrative
    parameters in millions of cubic metres and days.

    Args:
        river_inflow:
            The mean river inflow to the reservoir, per day.

        seasonality:
            The relative amplitude of the yearly inflow cycle.

        reservoir_capacity:
            The capacity of the reservoir.

        potential_et:
            The crop water requirement, the potential evapotranspiration of the irrigated area per day.

        canal_loss:
            The fraction of the release lost from the canals to the aquifer.

        pump_capacity:
            The largest groundwater pumping per day.

    Returns:
        system:
            The model, with the water stocks `scheme.reservoir`, `scheme.soil`, and `scheme.groundwater`, and the
            flows `scheme.inflow`, `scheme.evaporation`, `scheme.spill`, `scheme.irrigation`, `scheme.seepage`,
            `scheme.groundwater_pumping`, `scheme.crop_water`, `scheme.percolation`, and `scheme.baseflow` per day,
            over days. The observables of the model are `scheme.demand`, `scheme.release`, and `scheme.pumping`.

    Examples:
        >>> from skcomplex.models.engineering.processes import irrigation
        >>> irrigation(pump_capacity=0.0).elements[0].pump_capacity
        0.0
    """
    scheme = IrrigationSchemeElement(
        'scheme',
        river_inflow=river_inflow,
        seasonality=seasonality,
        reservoir_capacity=reservoir_capacity,
        potential_et=potential_et,
        canal_loss=canal_loss,
        pump_capacity=pump_capacity,
    )
    return System(elements=[scheme])

Components

class IrrigationSchemeElement(Element):
    """Irrigation scheme that moves water between a reservoir, the root zone of crops, and an aquifer."""

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        reservoir: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo), quantity='water')] = 50.0,
        soil: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo), quantity='water')] = 12.0,
        groundwater: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo), quantity='water')] = 100.0,
        river_inflow: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
        seasonality: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.9,
        reservoir_capacity: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 80.0,
        evaporation: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.001,
        spill_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        release_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
        canal_loss: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
        potential_et: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.4,
        field_capacity: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 20.0,
        target: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.7,
        refill_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 5.0,
        percolation: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.005,
        pump_capacity: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
        pump_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
        baseflow_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 500.0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.reservoir = reservoir
        self.soil = soil
        self.groundwater = groundwater
        self.river_inflow = river_inflow
        self.seasonality = seasonality
        self.reservoir_capacity = reservoir_capacity
        self.evaporation = evaporation
        self.spill_time = spill_time
        self.release_time = release_time
        self.canal_loss = canal_loss
        self.potential_et = potential_et
        self.field_capacity = field_capacity
        self.target = target
        self.refill_time = refill_time
        self.percolation = percolation
        self.pump_capacity = pump_capacity
        self.pump_time = pump_time
        self.baseflow_time = baseflow_time

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the flows of water through the scheme."""
        return IrrigationSchemeElement._irrigation_laws(s)

    @staticmethod
    def _demand(s: SimpleNamespace) -> Expression:
        """Return the water demand of the crops, their evapotranspiration and the refill of the root zone."""
        return clip(s.potential_et + (s.target * s.field_capacity - s.soil) / s.refill_time, 0)

    @staticmethod
    def _reservoir_release(s: SimpleNamespace) -> Expression:
        """Return the release of the reservoir, the demand before the canal losses, limited by the reservoir."""
        return minimum(IrrigationSchemeElement._demand(s) / (1 - s.canal_loss), s.reservoir / s.release_time)

    @staticmethod
    def _pumping(s: SimpleNamespace) -> Expression:
        """Return the pumping of groundwater, the demand the release leaves, limited by the pumps and the aquifer."""
        unmet = IrrigationSchemeElement._demand(s) - (1 - s.canal_loss) * IrrigationSchemeElement._reservoir_release(s)
        return minimum(clip(unmet, 0), minimum(s.pump_capacity, s.groundwater / s.pump_time))

    @staticmethod
    def _irrigation_laws(s: SimpleNamespace) -> list[BaseLaw]:
        """Return the water balance of a reservoir, the root zone of the irrigated crops, and the aquifer below them."""
        release = IrrigationSchemeElement._reservoir_release(s)
        crop_water = s.potential_et * lookup(s.soil / s.field_capacity, WATER_STRESS)
        spill = clip(s.reservoir - s.reservoir_capacity, 0) / s.spill_time
        return [
            TransferLaw(
                Source('river'),
                s.reservoir,
                RiceMillElement._season(s.river_inflow, s.seasonality, s.t, DAYS_PER_YEAR),
                label='inflow',
            ),
            TransferLaw(s.reservoir, Sink('atmosphere'), s.evaporation * s.reservoir, label='evaporation'),
            TransferLaw(s.reservoir, Sink('downstream'), spill, label='spill'),
            TransferLaw(s.reservoir, s.soil, (1 - s.canal_loss) * release, label='irrigation'),
            TransferLaw(s.reservoir, s.groundwater, s.canal_loss * release, label='seepage'),
            TransferLaw(s.groundwater, s.soil, IrrigationSchemeElement._pumping(s), label='groundwater_pumping'),
            TransferLaw(s.soil, Sink('atmosphere'), crop_water, label='crop_water'),
            TransferLaw(s.soil, s.groundwater, s.percolation * s.soil, label='percolation'),
            TransferLaw(s.groundwater, Sink('downstream'), s.groundwater / s.baseflow_time, label='baseflow'),
        ]

Equations

scheme.reservoir: diff(reservoir_scheme(t), t) == (seasonality_scheme*sin(0.01721420632103996*t) + 1)*river_inflow_scheme - evaporation_scheme*reservoir_scheme(t) + (canal_loss_scheme - 1)*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) - canal_loss_scheme*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) - where(-reservoir_capacity_scheme + reservoir_scheme(t) < 0, 0, -reservoir_capacity_scheme + reservoir_scheme(t))/spill_time_scheme
scheme.soil: diff(soil_scheme(t), t) == -percolation_scheme*soil_scheme(t) - (canal_loss_scheme - 1)*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) - potential_et_scheme*where(soil_scheme(t)/field_capacity_scheme < 0.0, 0.0, where(soil_scheme(t)/field_capacity_scheme < 0.6, 1.6666666666666667*soil_scheme(t)/field_capacity_scheme, 1.0)) + where(where(groundwater_scheme(t)/pump_time_scheme < pump_capacity_scheme, groundwater_scheme(t)/pump_time_scheme, pump_capacity_scheme) < where((canal_loss_scheme - 1)*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) + where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme) < 0, 0, (canal_loss_scheme - 1)*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) + where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)), where(groundwater_scheme(t)/pump_time_scheme < pump_capacity_scheme, groundwater_scheme(t)/pump_time_scheme, pump_capacity_scheme), where((canal_loss_scheme - 1)*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) + where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme) < 0, 0, (canal_loss_scheme - 1)*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) + where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)))
scheme.groundwater: diff(groundwater_scheme(t), t) == percolation_scheme*soil_scheme(t) + canal_loss_scheme*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) - groundwater_scheme(t)/baseflow_time_scheme - where(where(groundwater_scheme(t)/pump_time_scheme < pump_capacity_scheme, groundwater_scheme(t)/pump_time_scheme, pump_capacity_scheme) < where((canal_loss_scheme - 1)*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) + where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme) < 0, 0, (canal_loss_scheme - 1)*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) + where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)), where(groundwater_scheme(t)/pump_time_scheme < pump_capacity_scheme, groundwater_scheme(t)/pump_time_scheme, pump_capacity_scheme), where((canal_loss_scheme - 1)*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) + where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme) < 0, 0, (canal_loss_scheme - 1)*where(reservoir_scheme(t)/release_time_scheme < -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1), reservoir_scheme(t)/release_time_scheme, -where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)/(canal_loss_scheme - 1)) + where(potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme < 0, 0, potential_et_scheme + (field_capacity_scheme*target_scheme - soil_scheme(t))/refill_time_scheme)))