Routing

The routing module implements hillslope and channel routing algorithms.

Overview

The module provides three main routing components:

• Hillslope routing: Accumulates lateral flow contributions from upslope cells following D8 flow directions
• Channel routing: Routes water through the river network
• Reservoir routing: Optionally simulates reservoir storage dynamics with time-varying regulation and stage-discharge relationships

Hillslope and channel routing functions use Numba JIT compilation for high-performance execution.

Functions

Route lateral flow ONE STEP from upslope cells to immediate downstream neighbors.

This function performs ONE-STEP routing matching MATLAB’s hill_route.m behavior. Each cell receives flow ONLY from its immediate upstream neighbors, NOT from all upslope cells. Water moves gradually cell-by-cell over multiple timesteps.

CRITICAL: This is NOT cumulative routing! To move water from headwaters to outlets requires calling this function once per timestep for many timesteps.

PERFORMANCE: Uses Numba JIT compilation for speedup over pure Python loops.

Parameters:

Name	Type	Description	Default
lateral_flow	ndarray	2D array of lateral flow from each cell [m³/s]. Shape: (nrows, ncols). NaN values indicate no-data cells.	required
flow_direction	ndarray	2D array of flow directions in MOBIDIC notation (1-8) [dimensionless]. Shape: (nrows, ncols). Flow directions are standardized to MOBIDIC convention during preprocessing, regardless of original format.	required

Returns:

Type	Description
ndarray	Upstream contribution array with same shape as lateral_flow [m³/s].
ndarray	Each cell contains ONLY flow from immediate upstream neighbors (NOT its own flow).

Notes

Flow direction coding (D8) - MOBIDIC convention from MATLAB stack8point.m: MOBIDIC notation (1-8): 7 6 5 8 · 4 1 2 3

Direction number -> (row_offset, col_offset) in matrix coordinates: 1: (-1, -1) northwest, 2: (-1, 0) north, 3: (-1, 1) northeast, 4: (0, 1) east, 5: (1, 1) southeast, 6: (1, 0) south, 7: (1, -1) southwest, 8: (0, -1) west

Special values: 0: Outlet cell (no downstream neighbor) -1: Basin outlet marker (set by preprocessor at cell with max flow accumulation)

Examples:

>>> # Simple 3x3 grid - all cells flow to center
>>> lateral_flow = np.ones((3, 3)) * 0.1  # 0.1 m³/s from each cell
>>> # Flow directions in MOBIDIC notation (all point toward center at [1,1])
>>> flow_direction = np.array([[5, 6, 7],
...                            [4, 0, 8],  # center [1,1] is outlet (0 = no flow out)
...                            [3, 2, 1]])  # all 8 neighbors drain to center
>>> upstream = hillslope_routing(lateral_flow, flow_direction)
>>> upstream[1, 1]  # Center receives flow from all 8 neighbors (one-step)
0.8  # 8 neighbors x 0.1 m³/s each = 0.8 m³/s (does NOT include center's own 0.1)

Source code in mobidic/core/routing.py

def hillslope_routing(
    lateral_flow: np.ndarray,
    flow_direction: np.ndarray,
) -> np.ndarray:
    """
    Route lateral flow ONE STEP from upslope cells to immediate downstream neighbors.

    This function performs ONE-STEP routing matching MATLAB's hill_route.m behavior.
    Each cell receives flow ONLY from its immediate upstream neighbors, NOT from all
    upslope cells. Water moves gradually cell-by-cell over multiple timesteps.

    CRITICAL: This is NOT cumulative routing! To move water from headwaters to outlets
    requires calling this function once per timestep for many timesteps.

    PERFORMANCE: Uses Numba JIT compilation for speedup over pure Python loops.

    Args:
        lateral_flow: 2D array of lateral flow from each cell [m³/s].
            Shape: (nrows, ncols). NaN values indicate no-data cells.
        flow_direction: 2D array of flow directions in MOBIDIC notation (1-8) [dimensionless].
            Shape: (nrows, ncols). Flow directions are standardized to MOBIDIC convention
            during preprocessing, regardless of original format.

    Returns:
        Upstream contribution array with same shape as lateral_flow [m³/s].
        Each cell contains ONLY flow from immediate upstream neighbors (NOT its own flow).

    Notes:
        Flow direction coding (D8) - MOBIDIC convention from MATLAB stack8point.m:
            MOBIDIC notation (1-8):
            7  6  5
            8  ·  4
            1  2  3

        Direction number -> (row_offset, col_offset) in matrix coordinates:
            1: (-1, -1) northwest,  2: (-1, 0) north,    3: (-1, 1) northeast,
            4: (0, 1) east,         5: (1, 1) southeast, 6: (1, 0) south,
            7: (1, -1) southwest,   8: (0, -1) west

        Special values:
            0: Outlet cell (no downstream neighbor)
            -1: Basin outlet marker (set by preprocessor at cell with max flow accumulation)

    Examples:
        >>> # Simple 3x3 grid - all cells flow to center
        >>> lateral_flow = np.ones((3, 3)) * 0.1  # 0.1 m³/s from each cell
        >>> # Flow directions in MOBIDIC notation (all point toward center at [1,1])
        >>> flow_direction = np.array([[5, 6, 7],
        ...                            [4, 0, 8],  # center [1,1] is outlet (0 = no flow out)
        ...                            [3, 2, 1]])  # all 8 neighbors drain to center
        >>> upstream = hillslope_routing(lateral_flow, flow_direction)
        >>> upstream[1, 1]  # Center receives flow from all 8 neighbors (one-step)
        0.8  # 8 neighbors x 0.1 m³/s each = 0.8 m³/s (does NOT include center's own 0.1)
    """

    # Validate inputs
    if lateral_flow.shape != flow_direction.shape:
        raise ValueError(
            f"lateral_flow shape {lateral_flow.shape} must match flow_direction shape {flow_direction.shape}"
        )

    nrows, ncols = lateral_flow.shape

    # Initialize upstream contribution to ZERO (matching MATLAB line 19: uf=0*fl)
    # MATLAB hill_route does NOT include cell's own flow - only upstream contributions
    upstream_contribution = np.zeros_like(lateral_flow)

    # Call Numba-compiled kernel for maximum performance
    # Flow direction is always in MOBIDIC format (1-8) after preprocessing
    _hillslope_routing_kernel(
        lateral_flow,
        flow_direction,
        upstream_contribution,
        nrows,
        ncols,
    )

    return upstream_contribution

Route water through river network using linear reservoir method.

The linear routing method models each reach as a simple reservoir with exponential recession. Water is routed from upstream reaches to downstream reaches following the network topology, with contributions from lateral inflows.

The routing equation for each reach is

Q_out(t+dt) = C3 * Q_out(t) + C4 * qL

Where

C3 = exp(-dt/K) [recession coefficient] C4 = 1 - C3 [lateral inflow coefficient] K = lag_time_s [s] (lag time used as storage coefficient) qL = lateral inflow + integrated upstream contributions [m³/s]

PERFORMANCE: Uses Numba JIT compilation for significant speedup over pure Python loops.

Parameters:

Name	Type	Description	Default
network	GeoDataFrame | dict	River network GeoDataFrame with columns: - mobidic_id: Internal reach ID (0-indexed) - upstream_1, upstream_2: Upstream reach IDs (NaN if none) - downstream: Downstream reach ID (NaN if outlet) - calc_order: Calculation order (lower values processed first) - lag_time_s: Lag time [s] (used as storage coefficient K) OR dictionary with pre-processed topology (from Simulation._preprocess_network_topology): - ‘upstream_1_idx’: numpy array of first upstream indices - ‘upstream_2_idx’: numpy array of second upstream indices - ‘n_upstream’: numpy array of upstream counts - ‘sorted_reach_idx’: numpy array of reach indices sorted by calc_order - ‘K’: numpy array of storage coefficients - ‘n_reaches’: number of reaches	required
discharge_initial	ndarray	Initial discharge for each reach [m³/s]. Shape: (n_reaches,). Indexed by DataFrame index (not mobidic_id).	required
lateral_inflow	ndarray	Lateral inflow to each reach during this time step [m³/s]. Shape: (n_reaches,). Indexed by DataFrame index (not mobidic_id).	required
dt	float	Time step duration [s].	required

Returns:

Type	Description
tuple[ndarray, dict]	Tuple containing: - discharge_final: Discharge at end of time step [m³/s], shape (n_reaches,) - routing_state: Dictionary with routing diagnostics: - ‘C3’: Recession coefficients for each reach - ‘C4’: Lateral inflow coefficients for each reach - ‘qL_total’: Total inflow (lateral + upstream) for each reach [m³/s]

Raises:

Type	Description
ValueError	If network is missing required columns
ValueError	If array shapes don’t match network size

Notes

• Reaches are processed in calc_order to ensure upstream reaches are computed before downstream reaches
• For upstream contributions, the function computes the mean integral of the exponential decay: ∫[C3^t dt] over the time step
• Negative discharges are clipped to zero with a warning

Examples:

>>> # Simple 2-reach network: reach 0 flows into reach 1
>>> network = gpd.GeoDataFrame({
...     'mobidic_id': [0, 1],
...     'upstream_1': [np.nan, 0],
...     'upstream_2': [np.nan, np.nan],
...     'downstream': [1, np.nan],
...     'calc_order': [0, 1],
...     'lag_time_s': [3600.0, 7200.0],  # 1 hour and 2 hours
...     'geometry': [...]
... })
>>> Q_init = np.array([10.0, 5.0])  # m³/s
>>> qL = np.array([2.0, 1.0])  # m³/s lateral inflow
>>> Q_final, state = linear_channel_routing(network, Q_init, qL, dt=900)

Source code in mobidic/core/routing.py

def linear_channel_routing(
    network: GeoDataFrame | dict,
    discharge_initial: np.ndarray,
    lateral_inflow: np.ndarray,
    dt: float,
) -> tuple[np.ndarray, dict]:
    """
    Route water through river network using linear reservoir method.

    The linear routing method models each reach as a simple reservoir with exponential
    recession. Water is routed from upstream reaches to downstream reaches following
    the network topology, with contributions from lateral inflows.

    The routing equation for each reach is:
        Q_out(t+dt) = C3 * Q_out(t) + C4 * qL

    Where:
        C3 = exp(-dt/K)  [recession coefficient]
        C4 = 1 - C3      [lateral inflow coefficient]
        K = lag_time_s [s] (lag time used as storage coefficient)
        qL = lateral inflow + integrated upstream contributions [m³/s]

    PERFORMANCE: Uses Numba JIT compilation for significant speedup over pure Python loops.

    Args:
        network: River network GeoDataFrame with columns:
            - mobidic_id: Internal reach ID (0-indexed)
            - upstream_1, upstream_2: Upstream reach IDs (NaN if none)
            - downstream: Downstream reach ID (NaN if outlet)
            - calc_order: Calculation order (lower values processed first)
            - lag_time_s: Lag time [s] (used as storage coefficient K)
            OR dictionary with pre-processed topology (from Simulation._preprocess_network_topology):
            - 'upstream_1_idx': numpy array of first upstream indices
            - 'upstream_2_idx': numpy array of second upstream indices
            - 'n_upstream': numpy array of upstream counts
            - 'sorted_reach_idx': numpy array of reach indices sorted by calc_order
            - 'K': numpy array of storage coefficients
            - 'n_reaches': number of reaches
        discharge_initial: Initial discharge for each reach [m³/s].
            Shape: (n_reaches,). Indexed by DataFrame index (not mobidic_id).
        lateral_inflow: Lateral inflow to each reach during this time step [m³/s].
            Shape: (n_reaches,). Indexed by DataFrame index (not mobidic_id).
        dt: Time step duration [s].

    Returns:
        Tuple containing:
            - discharge_final: Discharge at end of time step [m³/s], shape (n_reaches,)
            - routing_state: Dictionary with routing diagnostics:
                - 'C3': Recession coefficients for each reach
                - 'C4': Lateral inflow coefficients for each reach
                - 'qL_total': Total inflow (lateral + upstream) for each reach [m³/s]

    Raises:
        ValueError: If network is missing required columns
        ValueError: If array shapes don't match network size

    Notes:
        - Reaches are processed in calc_order to ensure upstream reaches are
          computed before downstream reaches
        - For upstream contributions, the function computes the mean integral
          of the exponential decay: ∫[C3^t dt] over the time step
        - Negative discharges are clipped to zero with a warning

    Examples:
        >>> # Simple 2-reach network: reach 0 flows into reach 1
        >>> network = gpd.GeoDataFrame({
        ...     'mobidic_id': [0, 1],
        ...     'upstream_1': [np.nan, 0],
        ...     'upstream_2': [np.nan, np.nan],
        ...     'downstream': [1, np.nan],
        ...     'calc_order': [0, 1],
        ...     'lag_time_s': [3600.0, 7200.0],  # 1 hour and 2 hours
        ...     'geometry': [...]
        ... })
        >>> Q_init = np.array([10.0, 5.0])  # m³/s
        >>> qL = np.array([2.0, 1.0])  # m³/s lateral inflow
        >>> Q_final, state = linear_channel_routing(network, Q_init, qL, dt=900)
    """
    # Check if network is pre-processed dictionary (from Simulation class)
    if isinstance(network, dict):
        # Use pre-processed topology (fast path)
        logger.debug(f"Starting linear channel routing for {network['n_reaches']} reaches, dt={dt}s")

        n_reaches = network["n_reaches"]
        upstream_1_idx = network["upstream_1_idx"]
        upstream_2_idx = network["upstream_2_idx"]
        n_upstream = network["n_upstream"]
        sorted_reach_idx = network["sorted_reach_idx"]
        K = network["K"]

    else:
        # Process GeoDataFrame (slower path, for backwards compatibility)
        logger.debug(f"Starting linear channel routing for {len(network)} reaches, dt={dt}s")

        # Validate inputs
        required_cols = ["mobidic_id", "upstream_1", "upstream_2", "downstream", "calc_order", "lag_time_s"]
        missing_cols = [col for col in required_cols if col not in network.columns]
        if missing_cols:
            raise ValueError(f"Network missing required columns: {missing_cols}")

        n_reaches = len(network)

        # Create mapping from mobidic_id to DataFrame index
        mobidic_id_to_idx = {int(network.at[idx, "mobidic_id"]): idx for idx in network.index}

        # Pre-extract topology to numpy arrays (Strategy 2)
        upstream_1_idx = np.array(
            [mobidic_id_to_idx.get(int(uid), -1) if pd.notna(uid) else -1 for uid in network["upstream_1"]],
            dtype=np.int32,
        )
        upstream_2_idx = np.array(
            [mobidic_id_to_idx.get(int(uid), -1) if pd.notna(uid) else -1 for uid in network["upstream_2"]],
            dtype=np.int32,
        )
        n_upstream = np.array(
            [
                (1 if pd.notna(network.at[idx, "upstream_1"]) else 0)
                + (1 if pd.notna(network.at[idx, "upstream_2"]) else 0)
                for idx in network.index
            ],
            dtype=np.int32,
        )

        # Get sorted reach indices
        sorted_reach_idx = network.sort_values("calc_order").index.values.astype(np.int32)

        # Extract K (lag time as storage coefficient)
        K = network["lag_time_s"].values

    # Validate array sizes
    if len(discharge_initial) != n_reaches:
        raise ValueError(f"discharge_initial length {len(discharge_initial)} must match number of reaches {n_reaches}")

    if len(lateral_inflow) != n_reaches:
        raise ValueError(f"lateral_inflow length {len(lateral_inflow)} must match number of reaches {n_reaches}")

    if dt <= 0:
        raise ValueError(f"Time step dt must be positive, got {dt}")

    # Initialize output arrays
    discharge_final = np.zeros(n_reaches, dtype=np.float64)
    qL_total = np.zeros(n_reaches, dtype=np.float64)

    # Calculate routing coefficients for all reaches using lag_time_s as K
    # Note: Division by zero (K=0) and NaN (K=NaN) are handled in the kernel
    # (line 231: if np.isnan(K[ki]) or K[ki] <= 0), so we suppress warnings here
    with np.errstate(divide="ignore", invalid="ignore"):
        C3 = np.exp(-dt / K)  # Recession coefficient
        C4 = 1 - C3  # Lateral inflow coefficient

    # Call Numba-compiled kernel for maximum performance
    _linear_routing_kernel(
        discharge_initial,
        lateral_inflow,
        sorted_reach_idx,
        upstream_1_idx,
        upstream_2_idx,
        n_upstream,
        K,
        C3,
        C4,
        discharge_final,
        qL_total,
    )

    # Check for negative discharges
    negative_mask = discharge_final < 0
    if np.any(negative_mask):
        n_negative = np.sum(negative_mask)
        min_discharge = np.min(discharge_final[negative_mask])
        logger.warning(
            f"Linear routing produced {n_negative} negative discharges "
            f"(min={min_discharge:.6f} m³/s). Clipping to zero."
        )
        discharge_final = np.maximum(discharge_final, 0.0)

    # Prepare routing state dictionary
    routing_state = {
        "C3": C3,
        "C4": C4,
        "qL_total": qL_total,
    }

    logger.debug(
        f"Linear channel routing completed. "
        f"Discharge range: [{discharge_final.min():.3f}, {discharge_final.max():.3f}] m³/s"
    )

    return discharge_final, routing_state

Route water through reservoirs.

This function performs reservoir water balance calculations including: - Volume update based on inflows and outflows - Stage calculation from volume - Discharge calculation from stage using regulation curves - Sub-stepping for numerical stability - Interaction with soil water balance in reservoir basins

Matching MATLAB: reservoir_routing_include.m (lines 1-225)

Parameters:

Name	Type	Description	Default
reservoirs_data	list	List of Reservoir objects from preprocessing	required
reservoir_states	list[ReservoirState]	List of current ReservoirState objects	required
reach_discharge	ndarray	Array of reach discharges (m3/s)	required
surface_runoff	ndarray	2D grid of surface runoff rates (m/s)	required
lateral_flow	ndarray	2D grid of lateral flow rates (m/s)	required
soil_wg	ndarray	2D grid of gravitational water content (m)	required
soil_wg0	ndarray	2D grid of gravitational water capacity (m)	required
current_time	datetime	Current simulation time	required
dt	float	Time step (s)	required
cell_area	float	Grid cell area (m2)	required
base_substeps	int	Base number of sub-steps (default: N_SUBSTEP_RESERVOIR from constants)	N_SUBSTEP_RESERVOIR

Returns:

Type	Description
tuple[list[ReservoirState], ndarray, ndarray, ndarray, ndarray]	Tuple of: - Updated reservoir states - Modified reach discharge array - Modified surface runoff grid - Modified lateral flow grid - Modified soil_wg grid

Source code in mobidic/core/reservoir.py

def reservoir_routing(
    reservoirs_data: list,
    reservoir_states: list[ReservoirState],
    reach_discharge: np.ndarray,
    surface_runoff: np.ndarray,
    lateral_flow: np.ndarray,
    soil_wg: np.ndarray,
    soil_wg0: np.ndarray,
    current_time: datetime,
    dt: float,
    cell_area: float,
    base_substeps: int = N_SUBSTEP_RESERVOIR,
) -> tuple[list[ReservoirState], np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Route water through reservoirs.

    This function performs reservoir water balance calculations including:
    - Volume update based on inflows and outflows
    - Stage calculation from volume
    - Discharge calculation from stage using regulation curves
    - Sub-stepping for numerical stability
    - Interaction with soil water balance in reservoir basins

    Matching MATLAB: reservoir_routing_include.m (lines 1-225)

    Args:
        reservoirs_data: List of Reservoir objects from preprocessing
        reservoir_states: List of current ReservoirState objects
        reach_discharge: Array of reach discharges (m3/s)
        surface_runoff: 2D grid of surface runoff rates (m/s)
        lateral_flow: 2D grid of lateral flow rates (m/s)
        soil_wg: 2D grid of gravitational water content (m)
        soil_wg0: 2D grid of gravitational water capacity (m)
        current_time: Current simulation time
        dt: Time step (s)
        cell_area: Grid cell area (m2)
        base_substeps: Base number of sub-steps (default: N_SUBSTEP_RESERVOIR from constants)

    Returns:
        Tuple of:
            - Updated reservoir states
            - Modified reach discharge array
            - Modified surface runoff grid
            - Modified lateral flow grid
            - Modified soil_wg grid
    """
    # Create copies to avoid modifying inputs
    reach_discharge = reach_discharge.copy()
    surface_runoff = surface_runoff.copy()
    lateral_flow = lateral_flow.copy()
    soil_wg = soil_wg.copy()

    # Convert current_time to pandas Timestamp for comparison
    current_timestamp = pd.Timestamp(current_time)

    # Loop through all reservoirs (matching MATLAB line 16)
    for i, reservoir in enumerate(reservoirs_data):
        # Check if reservoir is active at current time
        if reservoir.date_start is not None and current_timestamp < reservoir.date_start:
            # Reservoir not yet active
            reservoir_states[i].outflow = 0.0
            reservoir_states[i].inflow = 0.0
            reservoir_states[i].volume = 0.0
            reservoir_states[i].stage = 0.0
            continue

        # Find current regulation period (matching MATLAB lines 19-21)
        # period_times is a dict: {"000": "2020-01-01T00:00:00", "001": "2020-06-01T00:00:00", ...}
        period_times = reservoir.period_times
        if period_times is None or len(period_times) == 0:
            logger.warning(f"Reservoir {reservoir.id} has no regulation periods, skipping")
            reservoir_states[i].outflow = 0.0
            reservoir_states[i].inflow = 0.0
            continue

        # Convert period_times to list of (index, timestamp) tuples
        period_list = []
        for key, time_str in period_times.items():
            period_idx = int(key)
            period_time = pd.Timestamp(time_str)
            if period_time <= current_timestamp:
                period_list.append((period_idx, period_time))

        if len(period_list) == 0:
            # No active period yet
            reservoir_states[i].outflow = 0.0
            reservoir_states[i].inflow = 0.0
            continue

        # Get the most recent period (matching MATLAB: iT = iT(end))
        period_list.sort(key=lambda x: x[1])  # Sort by time
        iT = period_list[-1][0]  # Get the index of the most recent period

        # Get regulation curves for current period (matching MATLAB lines 27-38)
        period_key = f"{iT:03d}"
        lawH = np.array(reservoir.stage_discharge_h[period_key])
        lawQ = np.array(reservoir.stage_discharge_q[period_key])

        # Calculate number of sub-steps (matching MATLAB lines 23-48)
        nsbgo = _calculate_substeps(
            lawH,
            lawQ,
            reservoir.stage_storage_curve,
            dt,
            base_substeps,
        )

        if nsbgo > 1:
            logger.debug(f"Reservoir {reservoir.id}: using {nsbgo} sub-steps")

        # Get reservoir basin pixels (linear indices)
        ibac = reservoir.basin_pixels
        if ibac is None or len(ibac) == 0:
            logger.warning(f"Reservoir {reservoir.id} has no basin pixels, skipping")
            continue

        # Convert linear indices to 2D indices (Fortran order, matching MATLAB)
        ibac_row, ibac_col = np.unravel_index(ibac, surface_runoff.shape, order="F")

        # Initialize accumulated outflow (will be averaged over sub-steps)
        total_outflow = 0.0

        # Loop over sub-steps (matching MATLAB line 78)
        for ir in range(nsbgo):
            # Calculate stage from volume (matching MATLAB lines 80-105)
            # Using stage-storage curve interpolation
            stage = _interpolate_stage_from_volume(
                reservoir.stage_storage_curve,
                reservoir_states[i].volume,
            )
            reservoir_states[i].stage = stage

            # Calculate outflow from stage using regulation curve (matching MATLAB lines 107-120)
            discharge = _interpolate_discharge_from_stage(lawH, lawQ, stage)

            # Calculate inflow from upstream reaches (matching MATLAB line 132)
            inlet_reaches = reservoir.inlet_reaches
            if inlet_reaches is not None and len(inlet_reaches) > 0:
                inflow = np.sum(reach_discharge[inlet_reaches])
            else:
                inflow = 0.0

            reservoir_states[i].inflow = inflow

            # Calculate lateral inflows from basin cells (matching MATLAB line 134)
            # Sum of (pir + pid) over basin cells, converted to discharge (m3/s)
            pir_basin = surface_runoff[ibac_row, ibac_col]
            pid_basin = lateral_flow[ibac_row, ibac_col]
            lateral_inflow = np.nansum(pir_basin + pid_basin) * cell_area

            # Calculate soil water change in basin (matching MATLAB line 135)
            # DV1 = -sum(Wg0(ibac) - Wg(ibac)) * cellsize^2
            # wg_basin = soil_wg[ibac_row, ibac_col]
            # wg0_basin = soil_wg0[ibac_row, ibac_col]
            # soil_water_change = -np.nansum(wg0_basin - wg_basin) * cell_area
            soil_water_change = 0.0  # Disable soil water interaction for now

            # Update volume (matching MATLAB lines 134-136)
            # DV = (Qin - Qout + lateral_inflow) * dt/nsbgo
            DV = (inflow - discharge + lateral_inflow) * dt / nsbgo
            new_volume = reservoir_states[i].volume + DV + soil_water_change

            # Handle negative volumes (matching MATLAB lines 137-148)
            if new_volume < 0:
                # Try to reduce outflow to keep volume non-negative
                if new_volume + discharge * dt / nsbgo > 0:
                    # Reduce discharge proportionally (matching MATLAB line 139)
                    discharge = discharge + new_volume / (dt / nsbgo)
                    new_volume = 0.0
                    # Refill soil to capacity (matching MATLAB line 141)
                    # soil_wg[ibac_row, ibac_col] = wg0_basin
                else:
                    # Not enough water even with zero discharge
                    # Adjust volume without changing discharge (matching MATLAB line 143)
                    new_volume = new_volume + discharge * dt / nsbgo - soil_water_change
                    discharge = 0.0
            else:
                # Volume is positive, refill soil to capacity (matching MATLAB line 147)
                # soil_wg[ibac_row, ibac_col] = wg0_basin
                pass

            reservoir_states[i].volume = new_volume

            # Accumulate outflow for averaging (matching MATLAB lines 162-165)
            total_outflow += discharge

        # Average outflow over sub-steps
        reservoir_states[i].outflow = total_outflow / nsbgo

        # Check for invalid outflow (matching MATLAB lines 175-177)
        if not np.isfinite(reservoir_states[i].outflow) or reservoir_states[i].outflow < 0:
            logger.error(f"Invalid outflow for reservoir {reservoir.id}: {reservoir_states[i].outflow}, setting to 0")
            reservoir_states[i].outflow = 0.0

        # Zero out surface runoff and lateral flow in reservoir basin (matching MATLAB lines 178-179)
        surface_runoff[ibac_row, ibac_col] = 0.0
        lateral_flow[ibac_row, ibac_col] = 0.0

        # Zero out discharge of upstream reaches (matching MATLAB lines 180-184)
        if inlet_reaches is not None and len(inlet_reaches) > 0:
            reach_discharge[inlet_reaches] = 0.0

    return reservoir_states, reach_discharge, surface_runoff, lateral_flow, soil_wg

State variables for a single reservoir.

Attributes:

Name	Type	Description
volume	float	Current reservoir volume (m3)
stage	float	Current water stage/level (m)
inflow	float	Current inflow from upstream reaches (m3/s)
outflow	float	Current outflow/release (m3/s)
withdrawal	float	Current withdrawal for water use (m3/s) (not implemented yet)

Source code in mobidic/core/reservoir.py

@dataclass
class ReservoirState:
    """State variables for a single reservoir.

    Attributes:
        volume: Current reservoir volume (m3)
        stage: Current water stage/level (m)
        inflow: Current inflow from upstream reaches (m3/s)
        outflow: Current outflow/release (m3/s)
        withdrawal: Current withdrawal for water use (m3/s) (not implemented yet)
    """

    volume: float
    stage: float
    inflow: float = 0.0
    outflow: float = 0.0
    withdrawal: float = 0.0

Design features

Hillslope routing

• One-step routing: Each call routes water ONE STEP to immediate downstream neighbors
• Gradual propagation: Water moves cell-by-cell over multiple timesteps
• D8 flow direction: Supports MOBIDIC notation (1-8)
• NaN handling: Robust handling of no-data cells
• Outlet detection: Handles outlet cells (flow_dir = 0 or -1)
• Numba acceleration: JIT-compiled kernel for maximum performance

Channel routing

Currently implements linear reservoir routing with the following features:

• Linear reservoir model: Simple exponential decay for each reach
• Topological routing: Processes reaches in correct upstream→downstream order
• Binary tree topology: Supports 0, 1, or 2 upstream tributaries per reach
• Mean integral method: Computes mean upstream discharge over time step
• Mass conservative: Total water balance preserved
• Configurable storage: Uses storage coefficient (K) calculated by the preprocessor from network attributes

Reservoir routing

• Volume-stage-discharge model: Simulates reservoir storage dynamics with stage-storage and stage-discharge relationships
• Cubic spline interpolation: Smooth stage-volume curves using cubic splines
• Time-varying regulation: Supports multiple regulation periods (e.g., seasonal winter/summer operations)
• Adaptive sub-stepping: Automatically divides timestep for numerical stability based on discharge variability
• Negative volume handling: Gracefully handles edge cases when reservoir would empty
• Basin zeroing: Automatically zeros surface runoff and lateral flow in reservoir basin pixels
• Network integration: Identifies inlet/outlet reaches, adds outflow to outlet lateral inflow, zeros inlet discharge

Routing equations

The linear reservoir routing equation for each reach is:

\[ \frac{dQ}{dt} = A(q_L + U \cdot Q - Q) \]

where \(Q\) is discharge, \(q_L\) is lateral inflow (surface + groundwater), \(A\) is a diagonal matrix of inverse characteristic times (\(1/K\)), and \(U\) is a binary topology matrix indicating tributary connections.

Linear reservoir model

For each reach:

\[ C_3 = \exp\left(-\frac{\Delta t}{K}\right) \quad \text{[recession coefficient]} \]

\[ C_4 = 1 - C_3 \quad \text{[lateral inflow coefficient]} \]

\[ Q_{\text{out}}(t + \Delta t) = C_3 \cdot Q_{\text{out}}(t) + C_4 \cdot q_L \]

Where:

• \(K\) = storage coefficient (lag time) [s]
• \(q_L\) = lateral inflow + integrated upstream contributions [m³/s]
• \(\Delta t\) = time step [s]

The mean integral of upstream discharge is computed as:

\[ \overline{Q} = \frac{q_L}{C_4} + \frac{q_L - Q_{\text{initial}} \cdot C_4}{\ln(C_3)} \]

Special cases:

• \(C_3 = 1\) (\(K \to \infty\)): No attenuation, \(\overline{Q} = q_L / C_4\)
• \(C_3 \approx 0\) (\(K \to 0\)): Instant routing, \(\overline{Q} = q_L / C_4\)

Reservoir routing model

For each reservoir at each timestep:

1. Volume update:

\[ V(t + \Delta t) = V(t) + (Q_{\text{in}} - Q_{\text{out}} - W) \cdot \Delta t \]

where \(V\) is volume [m³], \(Q_{\text{in}}\) is total inflow (upstream discharge + basin contributions) [m³/s], \(Q_{\text{out}}\) is regulated discharge [m³/s], and \(W\) is withdrawal [m³/s].

2. Stage calculation:

\[ h = f_{\text{stage}}(V) \]

Using cubic spline interpolation of the stage-storage curve.

3. Regulation period determination:

Based on the current date and regulation schedule, select the active regulation curve (e.g., “winter” or “summer”).

4. Discharge calculation:

\[ Q_{\text{out}} = f_{\text{discharge}}(h, \text{period}) \]

Using linear interpolation of the stage-discharge curve for the active period.

5. Sub-stepping:

If discharge variability is high, divide \(\Delta t\) into \(N\) sub-steps where \(N\) is determined by:

\[ N = 2^k \quad \text{where } k = \max \left(0, \log_2(\text{max discharge change} / \text{threshold})\right) \]

This ensures numerical stability by limiting the rate of volume change per sub-step.

6. Negative volume handling:

If \(V(t + \Delta t) < 0\), reduce \(Q_{\text{out}}\) such that \(V(t + \Delta t) = 0\):

\[ Q_{\text{out, adjusted}} = \frac{V(t) + Q_{\text{in}} \cdot \Delta t - W \cdot \Delta t}{\Delta t} \]