Simulation

The simulation module implements the main time-stepping loop of the MOBIDIC hydrological model, orchestrating water balance calculations, routing, and I/O operations.

Overview

The simulation engine coordinates:

• Input data loading: GIS preprocessing and meteorological forcing (station-based or raster-based)
• Meteorological forcing: Automatically detects station data (with spatial interpolation) or pre-interpolated raster data
• State initialization: Initial conditions for soil, surface, and channel states (supports warm start)
• Time-stepping loop: Sequential water balance and routing calculations
• Station-based interpolation: Grid interpolation using IDW or nearest neighbor with pre-computed weights
• Raster-based forcing: Direct sampling from pre-interpolated grids with grid alignment validation
• Interpolated meteo output: Optional export of interpolated grids for subsequent raster-based runs
• PET calculation: Simple 1 mm/day constant rate (energy balance not yet implemented)
• Results storage: Time series collection and state snapshots with automatic file chunking
• Output generation: NetCDF states (with chunking) and Parquet/CSV reports
• Restart capability: Load and resume from previously saved states

Current implementation: Includes soil water balance, routing (hillslope, channel, reservoir), and state/report I/O. Energy balance and groundwater models not yet implemented.

Classes

MOBIDIC simulation engine.

This class orchestrates the hydrological simulation, including: - Loading input data (GIS, meteorology) - Initializing state variables - Running the main time-stepping loop - Saving results

Examples:

>>> from mobidic import load_config, load_gisdata, Simulation, MeteoData
>>> config = load_config("Arno.yaml")
>>> gisdata = load_gisdata("Arno_gisdata.nc", "Arno_network.parquet")
>>> forcing = MeteoData.from_netcdf("Arno_meteo.nc")
>>> sim = Simulation(gisdata, forcing, config)
>>> results = sim.run("2020-01-01", "2020-12-31")
>>> results.save_report("Arno_discharge.parquet")

Source code in mobidic/core/simulation.py

class Simulation:
    """MOBIDIC simulation engine.

    This class orchestrates the hydrological simulation, including:
    - Loading input data (GIS, meteorology)
    - Initializing state variables
    - Running the main time-stepping loop
    - Saving results

    Examples:
        >>> from mobidic import load_config, load_gisdata, Simulation, MeteoData
        >>> config = load_config("Arno.yaml")
        >>> gisdata = load_gisdata("Arno_gisdata.nc", "Arno_network.parquet")
        >>> forcing = MeteoData.from_netcdf("Arno_meteo.nc")
        >>> sim = Simulation(gisdata, forcing, config)
        >>> results = sim.run("2020-01-01", "2020-12-31")
        >>> results.save_report("Arno_discharge.parquet")
    """

    def __init__(
        self,
        gisdata: Any,  # GISData object
        forcing: MeteoData | MeteoRaster,
        config: MOBIDICConfig,
    ):
        """Initialize simulation.

        Args:
            gisdata: Preprocessed GIS data (from load_gisdata or run_preprocessing)
            forcing: Meteorological forcing data as MeteoData (stations) or MeteoRaster (gridded)
            config: MOBIDIC configuration
        """
        self.gisdata = gisdata
        self.forcing = forcing
        self.config = config

        # Detect forcing type and set up method
        if isinstance(forcing, MeteoRaster):
            self.forcing_mode = "raster"
            logger.info("Using raster meteorological forcing (no interpolation)")
            # Validate grid alignment
            forcing.validate_grid_alignment(gisdata.metadata)
            # Set forcing extraction method
            self._get_forcing_fn = self._get_raster_forcing
            # No interpolation weights needed
            self._interpolation_weights = None
        else:  # MeteoData
            self.forcing_mode = "station"
            logger.info("Using station meteorological forcing (with spatial interpolation)")
            # Set forcing extraction method
            self._get_forcing_fn = self._interpolate_forcing
            # Interpolation weights will be precomputed later

        # Extract grid metadata
        self.nrows, self.ncols = gisdata.metadata["shape"]
        self.resolution = gisdata.metadata["resolution"]
        self.xllcorner = gisdata.metadata["xllcorner"]
        self.yllcorner = gisdata.metadata["yllcorner"]

        # Flow direction type from config
        self.flow_dir_type = config.raster_settings.flow_dir_type

        # Extract grids from gisdata
        self.dtm = gisdata.grids["dtm"]
        self.flow_dir = gisdata.grids["flow_dir"]
        self.flow_acc = gisdata.grids["flow_acc"]
        self.wc0 = gisdata.grids["Wc0"]
        self.wg0 = gisdata.grids["Wg0"]
        self.ks = gisdata.grids["ks"]
        self.alpsur = gisdata.grids["alpsur"]
        self.hillslope_reach_map = gisdata.hillslope_reach_map

        # Preprocess Wg0, Wc0, Wp0, and Ws0
        self.wc0 = self.wc0 * self.config.parameters.multipliers.Wc_factor
        self.wg0 = self.wg0 * self.config.parameters.multipliers.Wg_factor
        self.wp0 = np.zeros((self.nrows, self.ncols))
        self.ws0 = np.zeros((self.nrows, self.ncols))

        # Apply minumum limits
        self.wc0 = np.maximum(self.wc0, const.W_MIN)
        self.wg0 = np.maximum(self.wg0, const.W_MIN)

        # Transition factor between gravitational and capillary storage
        Wg_Wc_tr = self.config.parameters.multipliers.Wg_Wc_tr
        if Wg_Wc_tr >= 0:
            wtot = self.wc0 + self.wg0
            self.wg0 = np.minimum(Wg_Wc_tr * self.wg0, wtot)
            self.wc0 = wtot - self.wg0

        # Optional grids
        # These may be None if not provided: use get() to avoid KeyError
        self.kf = gisdata.grids.get("kf")
        self.gamma = gisdata.grids.get("gamma")
        self.kappa = gisdata.grids.get("kappa")
        self.beta = gisdata.grids.get("beta")
        self.alpha = gisdata.grids.get("alpha")

        # River network
        self.network = gisdata.network

        # Reservoirs (optional)
        self.reservoirs = getattr(gisdata, "reservoirs", None)
        if self.reservoirs is not None:
            logger.info(f"Reservoirs loaded: {len(self.reservoirs)} reservoirs")
        else:
            logger.debug("No reservoirs in gisdata")

        # Time step
        self.dt = config.simulation.timestep

        # Prepare parameter grids
        self.param_grids = self._prepare_grids()

        # Preprocess and cache network topology for fast routing
        self._network_topology = self._preprocess_network_topology()

        # Pre-compute interpolation weights for meteorological forcing (station mode only)
        # Currently, only precipitation is used in the main loop
        if self.forcing_mode == "station":
            self._interpolation_weights = self._precompute_interpolation_weights(["precipitation"])
        # else: self._interpolation_weights already set to None above

        # Time indices cache (will be populated in run() when simulation period is known)
        self._time_indices_cache = None

        # Initialize state
        self.state = None

        logger.info(
            f"Simulation initialized: grid={self.nrows}x{self.ncols}, "
            f"dt={self.dt}s, network={len(self.network)} reaches"
        )

    def _initial_state(self) -> SimulationState:
        """Initialize the initial simulation state variables with initial conditions from configuration.

        Returns:
            Initial simulation state
        """
        logger.info("Initial simulation state")

        # Initialize soil water content from configuration
        wcsat = self.config.initial_conditions.Wcsat  # Fraction of capillary saturation
        wgsat = self.config.initial_conditions.Wgsat  # Fraction of gravitational saturation
        ws_init = self.config.initial_conditions.Ws  # Initial surface water depth [m]

        # Initial values
        wc = self.wc0.copy() * wcsat
        wg = self.wg0.copy() * wgsat
        ws = np.full((self.nrows, self.ncols), ws_init)
        wp = self.wp0.copy()  # Initial plant/canopy water content [m]

        # Set NaN outside domain (lines 298-300 in mobidic_sid.m)
        wc[np.isnan(self.flow_acc)] = np.nan
        wg[np.isnan(self.flow_acc)] = np.nan
        ws[np.isnan(self.flow_acc)] = np.nan
        wp[np.isnan(self.flow_acc)] = np.nan

        # Initialize river discharge and lateral inflow
        discharge = np.zeros(len(self.network))
        lateral_inflow = np.zeros(len(self.network))

        # Initialize reservoir states if reservoirs exist
        reservoir_states = None
        if self.reservoirs is not None:
            from mobidic.core.reservoir import ReservoirState

            reservoir_states = []
            for reservoir in self.reservoirs.reservoirs:
                # Initialize with configured initial volume
                state = ReservoirState(
                    volume=reservoir.initial_volume,
                    stage=0.0,  # Will be calculated in first timestep
                    inflow=0.0,
                    outflow=0.0,
                )
                reservoir_states.append(state)
            logger.info(f"Initialized {len(reservoir_states)} reservoir states")

        logger.success(
            f"State initialized. Initial conditions (average): "
            f"Wc={np.nanmean(wc) * 1000:.1f} mm, "
            f"Wg={np.nanmean(wg) * 1000:.1f} mm, "
            f"Ws={np.nanmean(ws) * 1000:.1f} mm, "
            f"Wp={np.nanmean(wp) * 1000:.1f} mm"
        )

        return SimulationState(wc, wg, wp, ws, discharge, lateral_inflow, reservoir_states)

    def set_initial_state(
        self,
        state: SimulationState | None = None,
        state_file: str | Path | None = None,
        time_index: int = -1,
    ) -> None:
        """Set the initial simulation state from a previous simulation.

        This method allows restarting a simulation from a previously saved state,
        enabling warm starts, multi-stage simulations, or continuation after interruption.

        Args:
            state: SimulationState object to use as initial state. If provided, state_file is ignored.
            state_file: Path to NetCDF state file to load. Used if state is None.
            time_index: Time index to load from state file (default: -1 for last timestep).
                Only used when loading from state_file.

        Raises:
            ValueError: If neither state nor state_file is provided, or if state_file doesn't exist.

        Examples:
            >>> # Method 1: Set state directly from SimulationState object
            >>> sim.set_initial_state(state=previous_state)
            >>>
            >>> # Method 2: Load from state file (last timestep)
            >>> sim.set_initial_state(state_file="states.nc")
            >>>
            >>> # Method 3: Load specific timestep from state file
            >>> sim.set_initial_state(state_file="states.nc", time_index=10)
        """
        if state is not None:
            # Use provided state directly
            logger.info("Setting initial state from provided SimulationState object")
            self.state = state
            logger.success(
                f"Initial state set. State contains: "
                f"Wc={np.nanmean(state.wc) * 1000:.1f} mm, "
                f"Wg={np.nanmean(state.wg) * 1000:.1f} mm, "
                f"Ws={np.nanmean(state.ws) * 1000:.1f} mm, "
                f"Q_mean={np.mean(state.discharge):.3f} m³/s"
            )
        elif state_file is not None:
            # Load from state file
            state_path = Path(state_file)
            if not state_path.exists():
                raise ValueError(f"State file not found: {state_path}")

            logger.info(f"Loading initial state from file: {state_path}")

            from mobidic.io import load_state

            loaded_state, state_time, metadata = load_state(
                input_path=state_path,
                network_size=len(self.network),
                time_index=time_index,
                config=self.config,
                gisdata=self.gisdata,
            )

            self.state = loaded_state
            logger.success(
                f"Initial state loaded from {state_path} at time {state_time}. "
                f"State contains: "
                f"Wc={np.nanmean(loaded_state.wc) * 1000:.1f} mm, "
                f"Wg={np.nanmean(loaded_state.wg) * 1000:.1f} mm, "
                f"Ws={np.nanmean(loaded_state.ws) * 1000:.1f} mm, "
                f"Q_mean={np.mean(loaded_state.discharge):.3f} m³/s"
            )
        else:
            raise ValueError("Either 'state' or 'state_file' must be provided")

    def _interpolate_forcing(
        self,
        time: datetime,
        variable: str,
        weights_cache: dict[str, np.ndarray | None] | None = None,
        time_step_index: int | None = None,
    ) -> np.ndarray:
        """Interpolate meteorological forcing data to grid.

        Args:
            time: Current simulation time
            variable: Variable name ('precipitation', 'temperature_min', 'temperature_max', etc.)
            weights_cache: Optional pre-computed interpolation weights for performance.
                If None, weights will be computed on-the-fly.
            time_step_index: Optional time step index for accessing cached time indices.
                If None, time indices will be computed on-the-fly.

        Returns:
            2D grid of interpolated values
        """
        # Get station data for this variable
        if variable not in self.forcing.stations:
            raise ValueError(f"Variable '{variable}' not found in forcing data")

        var_stations = self.forcing.stations[variable]
        n_stations = len(var_stations)

        if n_stations == 0:
            raise ValueError(f"No stations found for variable '{variable}'")

        # Pre-allocate arrays (faster than list append)
        station_x = np.zeros(n_stations, dtype=np.float64)
        station_y = np.zeros(n_stations, dtype=np.float64)
        station_elevation = np.zeros(n_stations, dtype=np.float64)
        station_values = np.zeros(n_stations, dtype=np.float64)
        valid_mask = np.ones(n_stations, dtype=bool)

        # Use cached time indices if available
        if (
            self._time_indices_cache is not None
            and time_step_index is not None
            and variable in self._time_indices_cache
        ):
            time_indices_var = self._time_indices_cache[variable]

            if time_indices_var is not None:
                # Fast path: use pre-computed indices
                indices = time_indices_var[:, time_step_index]

                for i, station in enumerate(var_stations):
                    if len(station["time"]) == 0:
                        valid_mask[i] = False
                        continue

                    idx = indices[i]
                    station_x[i] = station["x"]
                    station_y[i] = station["y"]
                    station_elevation[i] = station["elevation"]
                    station_values[i] = station["data"][idx]
            else:
                # Fallback for variables with no cached indices (shouldn't happen)
                valid_mask[:] = False
        else:
            # Fallback: compute time indices on-the-fly (slower, used when cache is not available)
            target_time = pd.Timestamp(time)

            for i, station in enumerate(var_stations):
                if len(station["time"]) == 0:
                    valid_mask[i] = False
                    continue

                time_index = station["time"].get_indexer([target_time], method="nearest")[0]

                if time_index >= 0 and time_index < len(station["data"]):
                    station_x[i] = station["x"]
                    station_y[i] = station["y"]
                    station_elevation[i] = station["elevation"]
                    station_values[i] = station["data"][time_index]
                else:
                    valid_mask[i] = False

        # Filter out invalid stations
        if not valid_mask.any():
            raise ValueError(f"No valid data found for {variable} at time {time}")

        station_x = station_x[valid_mask]
        station_y = station_y[valid_mask]
        station_elevation = station_elevation[valid_mask]
        station_values = station_values[valid_mask]

        # Use single resolution value (assume square cells)
        resolution = self.resolution[0]

        # Choose interpolation method based on variable
        if variable == "precipitation":
            # Check if there's any rainfall before interpolating (MATLAB line 1241: if nansum(pp))
            # This optimization skips interpolation when all precipitation values are zero or NaN
            if np.nansum(station_values) == 0 or np.all(np.isnan(station_values)):
                # Return zero grid (matching MATLAB line 1256: p_i = Mzeros)
                grid_values = np.zeros_like(self.dtm)
            else:
                # Get precipitation interpolation method from config
                precip_interp = self.config.simulation.precipitation_interp

                if precip_interp == "Nearest":
                    # Use nearest neighbor for precipitation (MATLAB: pluviomap.m when thiessen=1)
                    grid_values = precipitation_interpolation(
                        station_x,
                        station_y,
                        station_values,
                        self.dtm,
                        self.xllcorner,
                        self.yllcorner,
                        resolution,
                    )
                else:  # "IDW"
                    # Use IDW without elevation correction for precipitation
                    # MATLAB mobidic_sid.m line 1246: uses tmww3 with expon=2
                    # tmww3 = tmww^3 where tmww = 1/dist^2, so tmww3 = 1/dist^6
                    # This gives stronger weight to nearby stations (power=6)
                    # Use cached weights if available
                    weights_matrix = weights_cache.get(variable) if weights_cache else None
                    grid_values = station_interpolation(
                        station_x,
                        station_y,
                        station_elevation,
                        station_values,
                        self.dtm,
                        self.xllcorner,
                        self.yllcorner,
                        resolution,
                        weights_matrix=weights_matrix,
                        apply_elevation_correction=False,  # switchregz=0 in MATLAB
                        power=6.0,  # tmww3 in MATLAB: (1/dist^2)^3 = 1/dist^6
                    )

                # Convert from mm/h to m/s
                # Station data is stored in mm/h (same as raster data)
                grid_values = grid_values / 1000.0 / 3600.0
        else:
            # Use IDW with elevation correction for temperature and other variables
            # MATLAB calc_forcing_day.m uses different settings per variable:
            # - Temperature (min/max): switchregz=1, expon=2 (elevation correction, power=2)
            # - Humidity: switchregz=0, expon=2 (no elevation correction, power=2)
            # - Wind: switchregz=0, expon=0.5 (no elevation correction, power=0.5)
            # - Radiation: switchregz=0, expon=2 (no elevation correction, power=2)
            # Currently only temperature is implemented with elevation correction
            # Use cached weights if available
            weights_matrix = weights_cache.get(variable) if weights_cache else None
            grid_values = station_interpolation(
                station_x,
                station_y,
                station_elevation,
                station_values,
                self.dtm,
                self.xllcorner,
                self.yllcorner,
                resolution,
                weights_matrix=weights_matrix,
                apply_elevation_correction=True,  # True for temperature, False for others
                power=2.0,  # 2.0 for most, 0.5 for wind
            )

        return grid_values

    def _get_raster_forcing(
        self,
        time: datetime,
        variable: str,
        weights_cache: dict[str, np.ndarray | None] | None = None,
        time_step_index: int | None = None,
    ) -> np.ndarray:
        """Extract meteorological forcing from raster data.

        Args:
            time: Current simulation time
            variable: Variable name ('precipitation', 'temperature', etc.)
            weights_cache: Unused (for signature compatibility with _interpolate_forcing)
            time_step_index: Unused (for signature compatibility with _interpolate_forcing)

        Returns:
            2D grid in simulation units (m/s for precip, degC for temperature)
        """
        # Get grid from raster (returns values in file units: mm/h for precip, degC for temp)
        grid = self.forcing.get_timestep(time, variable)

        # Convert units to match _interpolate_forcing output
        if variable == "precipitation":
            # Convert mm/h to m/s: divide by 1000 * 3600
            grid = grid / 1000.0 / 3600.0

        # Temperature and other variables are already in correct units (degC)
        return grid

    def _calculate_pet(self, time: datetime) -> np.ndarray:
        """Calculate potential evapotranspiration.

        Currently energy balance is not implemented -> uses MATLAB approach: constant 1 mm/day.
        MATLAB: etp = Mones/(1000*3600*24) [mobidic_sid.m line 332]

        Args:
            time: Current simulation time

        Returns:
            PET grid [m] over time step
        """
        # Use MATLAB approach: constant 1 mm/day when no energy balance
        pet = calculate_pet((self.nrows, self.ncols), self.dt, pet_rate_mm_day=1.0)

        return pet

    def _prepare_grids(self) -> dict[str, np.ndarray]:
        """Prepare grids for simulation.

        Returns:
            Dictionary of parameter grids
        """
        # Get parameters from config
        params = self.config.parameters

        # Create grids from scalar parameters or raster files
        # If raster exists, use it; otherwise use scalar value
        param_grids = {}

        # Rainfall fraction f0: fraction of time step without rain [-]
        # Time-dependent parameter from mobidic_sid.m line 213
        f0_value = const.F0_CONSTANT * (1 - np.exp(-self.dt / (24 * 3600) * np.log(const.F0_CONSTANT / 0.10)))
        param_grids["f0"] = np.full((self.nrows, self.ncols), f0_value)
        param_grids["f0"][np.isnan(self.dtm)] = np.nan

        # Hydraulic conductivity [m/s]
        ks_factor = self.config.parameters.multipliers.ks_factor
        if self.ks is not None:
            param_grids["ks"] = self.ks * ks_factor
        else:
            param_grids["ks"] = np.full((self.nrows, self.ncols), params.soil.ks) * ks_factor

        # Flow coefficients
        if self.gamma is not None:
            param_grids["gamma"] = self.gamma
        else:
            param_grids["gamma"] = np.full((self.nrows, self.ncols), params.soil.gamma)

        if self.kappa is not None:
            param_grids["kappa"] = self.kappa
        else:
            param_grids["kappa"] = np.full((self.nrows, self.ncols), params.soil.kappa)

        if self.beta is not None:
            param_grids["beta"] = self.beta
        else:
            param_grids["beta"] = np.full((self.nrows, self.ncols), params.soil.beta)

        if self.alpha is not None:
            param_grids["alpha"] = self.alpha
        else:
            param_grids["alpha"] = np.full((self.nrows, self.ncols), params.soil.alpha)

        # Channelized flow fraction cha [-]
        # From buildgis_mysql_include.m line 655-656
        param_grids["cha"] = self.flow_acc / np.nanmax(self.flow_acc)
        param_grids["cha"] = np.where(param_grids["cha"] > 0, param_grids["cha"], 0.0)

        # Surface alpha parameter alpsur
        param_grids["alpsur"] = self.alpsur * param_grids["alpha"]

        return param_grids

    def _precompute_interpolation_weights(self, variables: list[str] | None = None) -> dict[str, np.ndarray | None]:
        """Pre-compute IDW interpolation weights for each meteorological variable.

        Weights depend only on station geometry, not on values, so they can be
        computed once and reused for all timesteps. This provides significant
        performance improvement for long simulations.

        Args:
            variables: List of variable names to compute weights for.
                If None, computes weights for all variables in forcing data.
                Examples: ["precipitation"], ["precipitation", "temperature_min"]

        Returns:
            Dictionary mapping variable names to weight matrices (3D arrays).
            For variables using nearest neighbor (precipitation with "nearest" method),
            the value is None since weights aren't needed.
        """
        from mobidic.core.interpolation import compute_idw_weights

        # Determine which variables to compute weights for
        if variables is None:
            # Default: compute weights for all variables
            variables_to_process = list(self.forcing.stations.keys())
            logger.debug("Pre-computing interpolation weights for all meteorological variables")
        else:
            # Validate that requested variables exist in forcing data
            available_vars = set(self.forcing.stations.keys())
            invalid_vars = [v for v in variables if v not in available_vars]
            if invalid_vars:
                raise ValueError(
                    f"Variables {invalid_vars} not found in forcing data. Available variables: {sorted(available_vars)}"
                )
            variables_to_process = variables
            logger.debug(f"Pre-computing interpolation weights for selected variables: {variables_to_process}")

        weights_cache = {}
        resolution = self.resolution[0]  # Use single resolution value

        for variable in variables_to_process:
            var_stations = self.forcing.stations[variable]

            if len(var_stations) == 0:
                logger.debug(f"No stations for {variable}, skipping weight computation")
                weights_cache[variable] = None
                continue

            # Get station coordinates (assuming they don't change over time)
            station_x = np.array([s["x"] for s in var_stations])
            station_y = np.array([s["y"] for s in var_stations])

            # Determine interpolation settings based on variable type
            if variable == "precipitation":
                precip_interp = self.config.simulation.precipitation_interp
                if precip_interp == "Nearest":
                    # Nearest neighbor doesn't use weights matrix
                    logger.debug(f"{variable}: using nearest neighbor (no weights needed)")
                    weights_cache[variable] = None
                else:  # "IDW"
                    # IDW with power=6 for precipitation
                    logger.debug(f"{variable}: pre-computing IDW weights (power=6)")
                    weights_cache[variable] = compute_idw_weights(
                        station_x,
                        station_y,
                        self.dtm,
                        self.xllcorner,
                        self.yllcorner,
                        resolution,
                        power=6.0,
                    )
            elif variable in ["temperature_min", "temperature_max"]:
                # Temperature: power=2, elevation correction applied to values (not weights)
                logger.debug(f"{variable}: pre-computing IDW weights (power=2)")
                weights_cache[variable] = compute_idw_weights(
                    station_x,
                    station_y,
                    self.dtm,
                    self.xllcorner,
                    self.yllcorner,
                    resolution,
                    power=2.0,
                )
            elif variable == "wind_speed":
                # Wind: power=0.5
                logger.debug(f"{variable}: pre-computing IDW weights (power=0.5)")
                weights_cache[variable] = compute_idw_weights(
                    station_x,
                    station_y,
                    self.dtm,
                    self.xllcorner,
                    self.yllcorner,
                    resolution,
                    power=0.5,
                )
            else:
                # Other variables (humidity, radiation): power=2, no elevation correction
                logger.debug(f"{variable}: pre-computing IDW weights (power=2)")
                weights_cache[variable] = compute_idw_weights(
                    station_x,
                    station_y,
                    self.dtm,
                    self.xllcorner,
                    self.yllcorner,
                    resolution,
                    power=2.0,
                )

        logger.success(f"Interpolation weights cached for {len(weights_cache)} variables")
        return weights_cache

    def _precompute_time_indices(
        self, simulation_times: pd.DatetimeIndex, variables: list[str] | None = None
    ) -> dict[str, np.ndarray]:
        """Pre-compute time indices for all stations and all simulation timesteps.

        This method pre-computes the mapping from simulation timesteps to station data
        indices, avoiding repeated time lookups during the main simulation loop.
        This provides significant performance improvement for long simulations.

        Args:
            simulation_times: DatetimeIndex of all simulation timesteps
            variables: List of variable names to compute time indices for.
                If None, computes indices for all variables in forcing data.
                Examples: ["precipitation"], ["precipitation", "temperature_min"]

        Returns:
            Dictionary mapping variable names to (n_stations × n_timesteps) index arrays.
            For each variable and timestep, the array contains the index into the
            station's data array that is nearest to that timestep.
        """
        # Determine which variables to compute indices for
        if variables is None:
            # Default: compute indices for all variables
            variables_to_process = list(self.forcing.stations.keys())
            logger.info(f"Pre-computing time indices for {len(simulation_times)} timesteps (all variables)")
        else:
            # Validate that requested variables exist in forcing data
            available_vars = set(self.forcing.stations.keys())
            invalid_vars = [v for v in variables if v not in available_vars]
            if invalid_vars:
                raise ValueError(
                    f"Variables {invalid_vars} not found in forcing data. Available variables: {sorted(available_vars)}"
                )
            variables_to_process = variables
            logger.info(
                f"Pre-computing time indices for {len(simulation_times)} timesteps "
                f"({len(variables_to_process)} variables: {variables_to_process})"
            )

        time_indices = {}
        n_times = len(simulation_times)

        for variable in variables_to_process:
            var_stations = self.forcing.stations[variable]
            n_stations = len(var_stations)

            if n_stations == 0:
                logger.debug(f"No stations for {variable}, skipping time index computation")
                time_indices[variable] = None
                continue

            # Pre-allocate index array (n_stations × n_timesteps)
            indices = np.zeros((n_stations, n_times), dtype=np.int32)

            for i, station in enumerate(var_stations):
                if len(station["time"]) == 0:
                    # Station has no data - indices will be 0 but won't be used
                    continue

                # Convert station times to numpy array for searchsorted
                station_times = station["time"].values.astype("datetime64[ns]")
                sim_times = simulation_times.values.astype("datetime64[ns]")

                # Find insertion points (left side)
                left_idx = np.searchsorted(station_times, sim_times, side="left")

                # For each simulation time, find the nearest station time
                # We need to check both left_idx and left_idx-1 (if valid)
                max_idx = len(station_times) - 1

                for t in range(n_times):
                    left = left_idx[t]

                    # Clamp to valid range
                    if left > max_idx:
                        # Past the end - use last index
                        indices[i, t] = max_idx
                    elif left == 0:
                        # Before the start - use first index
                        indices[i, t] = 0
                    else:
                        # Check both neighbors and pick nearest
                        right = left
                        left_candidate = left - 1

                        left_dist = abs((sim_times[t] - station_times[left_candidate]).astype("timedelta64[ns]"))
                        right_dist = abs((sim_times[t] - station_times[right]).astype("timedelta64[ns]"))

                        if left_dist <= right_dist:
                            indices[i, t] = left_candidate
                        else:
                            indices[i, t] = right

            time_indices[variable] = indices
            logger.debug(f"Pre-computed {n_stations} x {n_times} time indices for {variable}")

        logger.debug(f"Time indices cached for {len(time_indices)} variables")
        return time_indices

    def _preprocess_network_topology(self) -> dict:
        """Preprocess network topology for fast routing.

        This method pre-extracts network topology to numpy arrays and caches them
        for use in every time step. This avoids repeated pandas operations during
        the main simulation loop.

        Returns:
            Dictionary containing pre-processed 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 (lag_time_s)
                - 'n_reaches': number of reaches
        """
        logger.debug("Preprocessing network topology for fast routing")

        network = self.network
        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
        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 by calc_order
        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

        topology = {
            "upstream_1_idx": upstream_1_idx,
            "upstream_2_idx": upstream_2_idx,
            "n_upstream": n_upstream,
            "sorted_reach_idx": sorted_reach_idx,
            "K": K,
            "n_reaches": n_reaches,
        }

        logger.success(f"Network topology preprocessed: {n_reaches} reaches cached")

        return topology

    def _accumulate_lateral_inflow(self, lateral_flow: np.ndarray) -> np.ndarray:
        """Accumulate lateral flow contributions to each reach.

        Following MATLAB's approach (mobidic_sid.m line 224):
        - Only accumulate from cells where ch > 0 (valid reach IDs)
        - Skip NaN cells (outside domain)
        - Skip -9999 cells (inside domain but cannot reach river network)

        Vectorized implementation using np.bincount for performance.

        Args:
            lateral_flow: 2D grid of lateral flow [m³/s]

        Returns:
            1D array of lateral inflow to each reach [m³/s]
        """
        n_reaches = len(self.network)

        # Flatten lateral flow and hillslope-reach mapping
        lateral_flow_flat = lateral_flow.ravel("F")
        hillslope_map_flat = self.hillslope_reach_map.ravel("F")

        # Create mask for valid cells: not NaN, >= 0, and finite lateral flow
        # MATLAB: ko = find(isfinite(zz) & (ch>0)); % contributing pixels
        valid_mask = np.isfinite(hillslope_map_flat) & (hillslope_map_flat >= 0) & np.isfinite(lateral_flow_flat)

        # Extract valid reach IDs and flow values
        valid_reach_ids = hillslope_map_flat[valid_mask].astype(np.int32)
        valid_flows = lateral_flow_flat[valid_mask]

        # Accumulate using bincount (vectorized, very fast)
        lateral_inflow = np.bincount(
            valid_reach_ids,
            weights=valid_flows,
            minlength=n_reaches,
        )

        return lateral_inflow

    def _should_save_state(self, step: int, current_time: datetime) -> bool:
        """Determine if state should be saved at current timestep.

        Args:
            step: Current timestep index (0-indexed)
            current_time: Current simulation datetime

        Returns:
            True if state should be saved, False otherwise
        """
        state_settings = self.config.output_states_settings

        if state_settings.output_states == "final":
            # Only save final state (handled after loop)
            return False
        elif state_settings.output_states == "all":
            # Save at intervals specified by output_interval
            if state_settings.output_interval is None:
                # No interval specified, save every timestep
                return True
            else:
                # Check if enough time has elapsed since start
                # We save at multiples of output_interval
                interval_steps = int(state_settings.output_interval / self.dt)
                return (step + 1) % interval_steps == 0
        elif state_settings.output_states == "list":
            # Save at specific datetimes specified in output_list
            # Convert datetime strings to datetime objects and check if current_time matches
            for date_str in state_settings.output_list:
                target_time = datetime.fromisoformat(date_str.replace(" ", "T"))
                if current_time == target_time:
                    return True
            return False
        else:
            return False

    def run(
        self,
        start_date: str | datetime,
        end_date: str | datetime,
        save_states_interval: int | None = None,
        save_report_interval: int | None = None,
    ) -> SimulationResults:
        """Run simulation.

        Args:
            start_date: Simulation start date (YYYY-MM-DD or datetime)
            end_date: Simulation end date (YYYY-MM-DD or datetime)
            save_states_interval: Interval for saving states [s]. If None, use config value.
            save_report_interval: Interval for saving reports [s]. If None, use config value.

        Returns:
            SimulationResults object containing time series and final state
        """

        logger.info("")
        logger.info("=" * 80)
        logger.info(f"MOBIDICpy v{__version__} - SIMULATION")
        logger.info("=" * 80)
        logger.info(f"Basin: {self.config.basin.id}")
        logger.info(f"Parameter set: {self.config.basin.paramset_id}")
        logger.info("")

        # Convert dates to datetime
        if isinstance(start_date, str):
            start_date = datetime.fromisoformat(start_date)
        if isinstance(end_date, str):
            end_date = datetime.fromisoformat(end_date)

        # Initialize state (skip if already set via set_initial_state())
        if self.state is None:
            self.state = self._initial_state()
        else:
            logger.info("Using pre-set initial state (skipping default initialization)")

        # Initialize results container
        results = SimulationResults(self.config, simulation=self)
        discharge_ts = []
        lateral_inflow_ts = []
        time_ts = []

        # Calculate number of time steps (inclusive of end_date)
        n_steps = int((end_date - start_date).total_seconds() / self.dt) + 1
        logger.info(f"Number of time steps: {n_steps}")

        # Pre-compute time indices for all simulation timesteps (station mode only)
        # Currently, only precipitation is used in the main loop
        simulation_times = pd.date_range(start=start_date, periods=n_steps, freq=f"{self.dt}s")
        if self.forcing_mode == "station":
            self._time_indices_cache = self._precompute_time_indices(simulation_times, variables=["precipitation"])
        else:
            self._time_indices_cache = None

        # Validate output_list if using list-based state output
        state_settings = self.config.output_states_settings
        if state_settings.output_states == "list" and state_settings.output_list:
            logger.debug("Validating output_list against simulation time steps")
            missing_states = []
            for date_str in state_settings.output_list:
                try:
                    target_time = datetime.fromisoformat(date_str.replace(" ", "T"))
                    # Check if target_time exists in simulation_times
                    if target_time not in simulation_times:
                        missing_states.append(date_str)
                except (ValueError, TypeError) as e:
                    logger.warning(f"Invalid datetime format in output_list: '{date_str}' - {e}")
                    missing_states.append(date_str)

            if missing_states:
                logger.warning(
                    f"The following {len(missing_states)} state(s) from output_list "
                    f"will NOT be saved (not found among simulation time steps):"
                )
                for missing_state in missing_states:
                    logger.warning(f"  - {missing_state}")
                logger.warning(
                    f"Simulation runs from {start_date.isoformat()} to {end_date.isoformat()} with timestep={self.dt}s"
                )

        # Get total soil capacity
        wtot0 = self.wc0 + self.wg0

        # Cell area for unit conversion [m²]
        cell_area = self.resolution[0] * self.resolution[1]

        # Create ko mask: contributing pixels only (matching MATLAB mobidic_sid.m:224)
        # ko = find(isfinite(zz) & (ch>0));  % contributing pixels
        ko = np.isfinite(self.dtm) & (self.hillslope_reach_map >= 0)
        ko = np.where(ko.ravel("F"))[0]
        logger.debug(f"Contributing pixels (ko): {len(ko)} of {self.nrows * self.ncols} total cells")

        # Initialize flow variables for hillslope routing feedback (matching MATLAB mobidic_sid.m:1621-1625)
        # flr_prev: surface runoff from previous timestep [m/s] - will be routed to create pir
        # fld_prev: lateral flow from previous timestep [m/s] - will be routed to create pid
        flr_prev = np.zeros((self.nrows, self.ncols))
        fld_prev = np.zeros((self.nrows, self.ncols))

        # Calculate progress logging interval: either 20 steps total or every 30 seconds
        # Use whichever results in fewer log messages
        steps_per_intervals = max(1, n_steps // 20)
        progress_step_interval = steps_per_intervals
        progress_time_interval = 30.0  # seconds

        # Prepare states directory and initialize state writer (if enabled)
        state_settings = self.config.output_states_settings
        state_writer = None

        # Only instantiate StateWriter if state output is enabled
        if state_settings.output_states not in [None, "None"]:
            states_dir = Path(self.config.paths.states)
            states_dir.mkdir(parents=True, exist_ok=True)

            reservoir_size = len(self.reservoirs) if self.reservoirs is not None else 0
            state_writer = StateWriter(
                output_path=states_dir / "states.nc",
                grid_metadata=self.gisdata.metadata,
                network_size=len(self.network),
                output_states=self.config.output_states,
                flushing=state_settings.flushing,
                max_file_size=state_settings.max_file_size,
                add_metadata={
                    "basin_id": self.config.basin.id,
                    "paramset_id": self.config.basin.paramset_id,
                },
                reservoir_size=reservoir_size,
            )
            logger.info(f"State output enabled: {state_settings.output_states}")
        else:
            logger.info("State output disabled (output_states=None)")

        # Initialize meteo writer (if enabled)
        meteo_writer = None
        if self.config.output_forcing_data.meteo_data:
            output_dir = Path(self.config.paths.output)
            output_dir.mkdir(parents=True, exist_ok=True)

            from mobidic.io import MeteoWriter

            # Define which variables to save
            meteo_variables = ["precipitation"]

            meteo_writer = MeteoWriter(
                output_path=output_dir / "meteo_forcing.nc",
                grid_metadata=self.gisdata.metadata,
                variables=meteo_variables,
                add_metadata={
                    "basin_id": self.config.basin.id,
                    "paramset_id": self.config.basin.paramset_id,
                    "start_date": start_date.isoformat(),
                    "end_date": end_date.isoformat(),
                },
            )
            logger.info(f"Meteo data output enabled: {meteo_variables}")
        else:
            logger.info("Interpolated meteo data output disabled")

        # Main time loop
        logger.info("Starting simulation main loop")
        current_time = start_date
        simulation_start_time = time.time()
        last_log_time = simulation_start_time
        for step in range(n_steps):
            logger.debug(f"Time step {step + 1}/{n_steps}: {current_time}")

            # 1. Get meteorological forcing (interpolation for stations, direct sampling for rasters)
            try:
                precip = self._get_forcing_fn(
                    current_time,
                    "precipitation",
                    weights_cache=self._interpolation_weights,
                    time_step_index=step,
                )
            except (KeyError, ValueError):
                logger.warning(f"Precipitation data not found for {current_time}, using zero")
                precip = np.zeros((self.nrows, self.ncols))

            # Log precipitation statistics (precip is in m/s, convert to mm/h for readability)
            precip_mm_h = precip * 1000.0 * 3600.0  # m/s -> mm/h
            precip_valid = precip_mm_h[np.isfinite(precip_mm_h)]
            if len(precip_valid) > 0:
                logger.debug(
                    f"Precipitation: mean={np.mean(precip_valid):.4f} mm/h, max={np.max(precip_valid):.4f} mm/h"
                )

            # 2. Calculate PET
            pet = self._calculate_pet(current_time)

            # 3. Save interpolated meteorological data (if enabled)
            if meteo_writer is not None:
                meteo_writer.append(current_time, precipitation=precip)

            # 4. Hillslope routing of previous timestep's flows (matching MATLAB mobidic_sid.m:1621-1625)
            # This must happen BEFORE soil mass balance to provide upstream contributions
            logger.debug("Routing previous timestep's flows through hillslope")

            # Convert to discharge for routing
            flr_prev_discharge = flr_prev * cell_area
            fld_prev_discharge = fld_prev * cell_area

            # Route through hillslope
            pir_discharge = hillslope_routing(flr_prev_discharge, self.flow_dir)
            pid_discharge = hillslope_routing(fld_prev_discharge, self.flow_dir)
            logger.debug("Hillslope routing completed")

            # Convert back to rates
            pir = pir_discharge / cell_area
            pid = pid_discharge / cell_area

            # 5. Soil water balance (cell-by-cell)
            logger.debug("Computing soil water balance")

            # Convert precipitation to depth over time step [m]
            precip_depth = precip * self.dt

            # Flatten 2D arrays to 1D (MATLAB: uses linear indexing with ko)
            # Extract only contributing pixels (ko) for processing (matching MATLAB mobidic_sid.m:733-736)
            wc_flat_full = self.state.wc.ravel("F")
            wg_flat_full = self.state.wg.ravel("F")
            wp_flat_full = self.state.wp.ravel("F") if self.state.wp is not None else None
            ws_flat_full = self.state.ws.ravel("F")

            # Extract ko cells only (matching MATLAB: Wc_ko = Wc(ko))
            wc_flat = wc_flat_full[ko]
            wg_flat = wg_flat_full[ko]
            wp_flat = wp_flat_full[ko] if wp_flat_full is not None else None
            ws_flat = ws_flat_full[ko]
            wc0_flat = self.wc0.ravel("F")[ko]
            wg0_flat = self.wg0.ravel("F")[ko]
            wtot0_flat = wtot0.ravel("F")[ko]
            precip_flat = precip_depth.ravel("F")[ko]
            pet_flat = pet.ravel("F")[ko] * self.dt
            # Multiply rate parameters by dt (matching MATLAB mobidic_sid.m:1681-1682)
            ks_flat = self.param_grids["ks"].ravel("F")[ko] * self.dt
            gamma_flat = self.param_grids["gamma"].ravel("F")[ko] * self.dt
            kappa_flat = self.param_grids["kappa"].ravel("F")[ko] * self.dt
            beta_flat = self.param_grids["beta"].ravel("F")[ko] * self.dt
            cha_flat = self.param_grids["cha"].ravel("F")[ko]
            f0_flat = self.param_grids["f0"].ravel("F")[ko]
            alpsur_flat = self.param_grids["alpsur"].ravel("F")[ko] * self.dt

            # Prepare routed flows from previous timestep (matching MATLAB mobidic_sid.m:1680)
            # Convert from [m/s] to [m] by multiplying by dt
            pir_flat = pir.ravel("F")[ko] * self.dt
            pid_flat = pid.ravel("F")[ko] * self.dt

            # Call soil_mass_balance with flattened arrays
            # Parameters have been pre-multiplied by dt above (lines 609-612, 616)
            (
                wc_out_flat,
                wg_out_flat,
                wp_out_flat,
                ws_out_flat,
                surface_runoff_flat,
                lateral_flow_depth_flat,
                et_flat,
                percolation_flat,
                capillary_flux_flat,
                wg_before_flat,
            ) = soil_mass_balance(
                wc=wc_flat,
                wc0=wc0_flat,
                wg=wg_flat,
                wg0=wg0_flat,
                wp=wp_flat,
                wp0=None,  # Plant reservoir disabled for now
                ws=ws_flat,
                ws0=None,
                wtot0=wtot0_flat,
                precipitation=precip_flat,
                surface_runoff_in=pir_flat,  # Routed surface runoff from previous timestep
                lateral_flow_in=pid_flat,  # Routed lateral flow from previous timestep
                potential_et=pet_flat,
                hydraulic_conductivity=ks_flat,
                hydraulic_conductivity_min=None,
                hydraulic_conductivity_max=None,
                channelized_fraction=cha_flat,
                surface_flow_exp=np.zeros_like(wc_flat),  # Not used
                lateral_flow_coeff=beta_flat,
                percolation_coeff=gamma_flat,
                absorption_coeff=kappa_flat,
                rainfall_fraction=f0_flat,
                et_shape=3.0,
                capillary_rise_enabled=False,
                test_mode=False,
                alpsur=alpsur_flat,
            )

            # Write ko results back to full grids (matching MATLAB: Wc(ko) = Wc_ko)
            wc_flat_full[ko] = wc_out_flat
            wg_flat_full[ko] = wg_out_flat
            if wp_out_flat is not None:
                wp_flat_full[ko] = wp_out_flat
            ws_flat_full[ko] = ws_out_flat

            # Initialize full output arrays for fluxes (preserve NaN for cells outside domain)
            surface_runoff_full = np.full(self.nrows * self.ncols, np.nan)
            lateral_flow_depth_full = np.full(self.nrows * self.ncols, np.nan)
            # Set non-ko cells within basin to zero
            surface_runoff_full[np.isfinite(self.dtm.ravel("F"))] = 0.0
            lateral_flow_depth_full[np.isfinite(self.dtm.ravel("F"))] = 0.0
            # Update ko cells with computed values
            surface_runoff_full[ko] = surface_runoff_flat
            lateral_flow_depth_full[ko] = lateral_flow_depth_flat

            # Reshape outputs back to 2D (Fortran order to match MATLAB)
            self.state.wc = wc_flat_full.reshape((self.nrows, self.ncols), order="F")
            self.state.wg = wg_flat_full.reshape((self.nrows, self.ncols), order="F")
            self.state.wp = (
                wp_flat_full.reshape((self.nrows, self.ncols), order="F") if wp_flat_full is not None else None
            )
            self.state.ws = ws_flat_full.reshape((self.nrows, self.ncols), order="F")
            surface_runoff = surface_runoff_full.reshape((self.nrows, self.ncols), order="F")
            lateral_flow_depth = lateral_flow_depth_full.reshape((self.nrows, self.ncols), order="F")

            # 5. Convert outputs from depth [m] to rate [m/s] (matching MATLAB mobidic_sid.m:1684)
            # flr: surface runoff rate [m/s] - analogous to MATLAB flr after division by dt
            # fld: lateral flow rate [m/s] - analogous to MATLAB fld after division by dt
            flr = surface_runoff / self.dt
            fld = lateral_flow_depth / self.dt

            # 6. Accumulate lateral inflow to reaches from surface runoff (matching MATLAB glob_route_day.m)
            # Convert to discharge for accumulation
            flr_discharge = flr * cell_area
            lateral_inflow = self._accumulate_lateral_inflow(flr_discharge)
            logger.debug("Lateral inflow accumulation completed")

            # Store lateral inflow in state
            self.state.lateral_inflow = lateral_inflow

            # Zero out flr for ALL cells that contributed to reaches (matching MATLAB glob_route_day.m line 33)
            # This prevents double-counting - flows from all contributing cells are consumed after accumulation
            flr[self.hillslope_reach_map >= 0] = 0.0
            # Note: fld is NOT zeroed - lateral flow continues to route between hillslope cells

            # Store flr and fld for next timestep's routing (at step 3)
            flr_prev = flr.copy()
            fld_prev = fld.copy()

            # 6. Reservoir routing (if reservoirs exist)
            # Prepare network topology for routing (may be modified if reservoirs exist)
            routing_network = self._network_topology

            if self.reservoirs is not None and self.state.reservoir_states is not None:
                logger.debug("Computing reservoir routing")

                # Call reservoir routing
                (
                    self.state.reservoir_states,
                    self.state.discharge,
                    flr,
                    fld,
                    self.state.wg,
                ) = reservoir_routing(
                    reservoirs_data=self.reservoirs.reservoirs,
                    reservoir_states=self.state.reservoir_states,
                    reach_discharge=self.state.discharge,
                    surface_runoff=flr,
                    lateral_flow=fld,
                    soil_wg=self.state.wg,
                    soil_wg0=self.wg0,
                    current_time=current_time,
                    dt=self.dt,
                    cell_area=cell_area,
                )

                # Add reservoir outflows to lateral inflow of outlet reaches (matching MATLAB glob_route_day.m:46-53)
                for i, reservoir in enumerate(self.reservoirs.reservoirs):
                    outlet_reach = reservoir.outlet_reach
                    if outlet_reach is not None:
                        lateral_inflow[outlet_reach] += self.state.reservoir_states[i].outflow

                # Clear upstream connections for outlet reaches (matching MATLAB glob_route_day.m:48-50)
                # This prevents double-counting: the reservoir has already collected upstream volumes
                # Create modified topology with cleared upstream connections for outlet reaches
                routing_network = self._network_topology.copy()
                routing_network["upstream_1_idx"] = self._network_topology["upstream_1_idx"].copy()
                routing_network["upstream_2_idx"] = self._network_topology["upstream_2_idx"].copy()
                routing_network["n_upstream"] = self._network_topology["n_upstream"].copy()

                for reservoir in self.reservoirs.reservoirs:
                    outlet_reach = reservoir.outlet_reach
                    if outlet_reach is not None:
                        # Clear upstream connections (MATLAB: ret(j).ramimonte=[nan nan])
                        routing_network["upstream_1_idx"][outlet_reach] = -1
                        routing_network["upstream_2_idx"][outlet_reach] = -1
                        # Set upstream count to 0 (MATLAB: ram_fin(j) = 0)
                        routing_network["n_upstream"][outlet_reach] = 0

            # 7. Channel routing
            logger.debug("Computing channel routing")
            self.state.discharge, routing_state = linear_channel_routing(
                network=routing_network,
                discharge_initial=self.state.discharge,
                lateral_inflow=lateral_inflow,
                dt=self.dt,
            )

            # 8. Store results
            discharge_ts.append(self.state.discharge.copy())
            lateral_inflow_ts.append(lateral_inflow.copy())
            time_ts.append(current_time)

            # 9. Save intermediate states if configured
            if state_writer is not None and self._should_save_state(step, current_time):
                logger.debug(f"Appending state to buffer at step {step + 1}/{n_steps}")
                state_writer.append_state(self.state, current_time)

            # Log progress based on step interval or time interval (whichever comes first)
            current_wall_time = time.time()
            time_since_last_log = current_wall_time - last_log_time
            is_step_interval = (step + 1) % progress_step_interval == 0
            is_time_interval = time_since_last_log >= progress_time_interval
            is_final_step = step == n_steps - 1

            if is_step_interval or is_time_interval or is_final_step:
                progress_bar = _create_progress_bar(step + 1, n_steps, bar_length=20)
                q_mean = np.mean(self.state.discharge)
                q_max = np.max(self.state.discharge)
                date_str = current_time.strftime("%Y-%m-%d %H:%M")

                logger.info(
                    f"{progress_bar} {step + 1}/{n_steps} | Simulation time: {date_str} | "
                    f"Q_mean={q_mean:.3f} m³/s | Q_max={q_max:.3f} m³/s"
                )
                last_log_time = current_wall_time

            # Advance time
            current_time += timedelta(seconds=self.dt)

        # Store results
        results.time_series["discharge"] = np.array(discharge_ts)
        results.time_series["lateral_inflow"] = np.array(lateral_inflow_ts)
        results.time_series["time"] = time_ts
        results.final_state = self.state

        # Calculate elapsed time
        simulation_end_time = time.time()
        elapsed_seconds = simulation_end_time - simulation_start_time
        elapsed_minutes = elapsed_seconds / 60
        elapsed_hours = elapsed_minutes / 60

        # Format elapsed time message
        if elapsed_hours >= 1:
            elapsed_str = f"{elapsed_hours:.2f} hours"
        elif elapsed_minutes >= 1:
            elapsed_str = f"{elapsed_minutes:.2f} minutes"
        else:
            elapsed_str = f"{elapsed_seconds:.2f} seconds"

        logger.success(f"Simulation completed: {n_steps} time steps in {elapsed_str}")

        # Auto-save reports based on configuration
        output_dir = Path(self.config.paths.output)
        output_dir.mkdir(parents=True, exist_ok=True)

        # Determine reach selection
        settings = self.config.output_report_settings
        reach_selection = settings.reach_selection
        selected_reaches = None
        reach_file = None

        if settings.reach_selection == "file":
            reach_file = settings.sel_file
        elif settings.reach_selection == "list":
            selected_reaches = settings.sel_list

        # Get output format from configuration
        output_format = self.config.output_report_settings.output_format
        file_extension = "csv" if output_format.lower() == "csv" else "parquet"

        # Save discharge report if enabled
        if self.config.output_report.discharge:
            start_str = start_date.strftime("%Y%m%d")
            end_str = end_date.strftime("%Y%m%d")
            discharge_path = output_dir / f"discharge_{start_str}_{end_str}.{file_extension}"
            logger.info("Exporting discharges")
            results.save_report(
                output_path=discharge_path,
                reach_selection=reach_selection,
                selected_reaches=selected_reaches,
                reach_file=reach_file,
                output_format=output_format,
            )

        # Save lateral inflow report if enabled
        if self.config.output_report.lateral_inflow:
            lateral_inflow_path = output_dir / f"lateral_inflow_{start_str}_{end_str}.{file_extension}"
            logger.info("Exporting lateral inflows")
            results.save_lateral_inflow_report(
                output_path=lateral_inflow_path,
                reach_selection=reach_selection,
                selected_reaches=selected_reaches,
                reach_file=reach_file,
                output_format=output_format,
            )

        # Save final state if enabled (for "final" mode, only save the last state)
        if state_writer is not None:
            if state_settings.output_states == "final":
                final_time = results.time_series["time"][-1]
                logger.info("Saving final state")
                state_writer.append_state(self.state, final_time)

            # Close the state writer (flushes any remaining buffered states)
            state_writer.close()

        # Close the meteo writer (writes all buffered data to NetCDF)
        if meteo_writer is not None:
            meteo_writer.close()

        return results

__init__(gisdata, forcing, config)

Initialize simulation.

Parameters:

Name	Type	Description	Default
gisdata	Any	Preprocessed GIS data (from load_gisdata or run_preprocessing)	required
forcing	MeteoData | MeteoRaster	Meteorological forcing data as MeteoData (stations) or MeteoRaster (gridded)	required
config	MOBIDICConfig	MOBIDIC configuration	required

Source code in mobidic/core/simulation.py

def __init__(
    self,
    gisdata: Any,  # GISData object
    forcing: MeteoData | MeteoRaster,
    config: MOBIDICConfig,
):
    """Initialize simulation.

    Args:
        gisdata: Preprocessed GIS data (from load_gisdata or run_preprocessing)
        forcing: Meteorological forcing data as MeteoData (stations) or MeteoRaster (gridded)
        config: MOBIDIC configuration
    """
    self.gisdata = gisdata
    self.forcing = forcing
    self.config = config

    # Detect forcing type and set up method
    if isinstance(forcing, MeteoRaster):
        self.forcing_mode = "raster"
        logger.info("Using raster meteorological forcing (no interpolation)")
        # Validate grid alignment
        forcing.validate_grid_alignment(gisdata.metadata)
        # Set forcing extraction method
        self._get_forcing_fn = self._get_raster_forcing
        # No interpolation weights needed
        self._interpolation_weights = None
    else:  # MeteoData
        self.forcing_mode = "station"
        logger.info("Using station meteorological forcing (with spatial interpolation)")
        # Set forcing extraction method
        self._get_forcing_fn = self._interpolate_forcing
        # Interpolation weights will be precomputed later

    # Extract grid metadata
    self.nrows, self.ncols = gisdata.metadata["shape"]
    self.resolution = gisdata.metadata["resolution"]
    self.xllcorner = gisdata.metadata["xllcorner"]
    self.yllcorner = gisdata.metadata["yllcorner"]

    # Flow direction type from config
    self.flow_dir_type = config.raster_settings.flow_dir_type

    # Extract grids from gisdata
    self.dtm = gisdata.grids["dtm"]
    self.flow_dir = gisdata.grids["flow_dir"]
    self.flow_acc = gisdata.grids["flow_acc"]
    self.wc0 = gisdata.grids["Wc0"]
    self.wg0 = gisdata.grids["Wg0"]
    self.ks = gisdata.grids["ks"]
    self.alpsur = gisdata.grids["alpsur"]
    self.hillslope_reach_map = gisdata.hillslope_reach_map

    # Preprocess Wg0, Wc0, Wp0, and Ws0
    self.wc0 = self.wc0 * self.config.parameters.multipliers.Wc_factor
    self.wg0 = self.wg0 * self.config.parameters.multipliers.Wg_factor
    self.wp0 = np.zeros((self.nrows, self.ncols))
    self.ws0 = np.zeros((self.nrows, self.ncols))

    # Apply minumum limits
    self.wc0 = np.maximum(self.wc0, const.W_MIN)
    self.wg0 = np.maximum(self.wg0, const.W_MIN)

    # Transition factor between gravitational and capillary storage
    Wg_Wc_tr = self.config.parameters.multipliers.Wg_Wc_tr
    if Wg_Wc_tr >= 0:
        wtot = self.wc0 + self.wg0
        self.wg0 = np.minimum(Wg_Wc_tr * self.wg0, wtot)
        self.wc0 = wtot - self.wg0

    # Optional grids
    # These may be None if not provided: use get() to avoid KeyError
    self.kf = gisdata.grids.get("kf")
    self.gamma = gisdata.grids.get("gamma")
    self.kappa = gisdata.grids.get("kappa")
    self.beta = gisdata.grids.get("beta")
    self.alpha = gisdata.grids.get("alpha")

    # River network
    self.network = gisdata.network

    # Reservoirs (optional)
    self.reservoirs = getattr(gisdata, "reservoirs", None)
    if self.reservoirs is not None:
        logger.info(f"Reservoirs loaded: {len(self.reservoirs)} reservoirs")
    else:
        logger.debug("No reservoirs in gisdata")

    # Time step
    self.dt = config.simulation.timestep

    # Prepare parameter grids
    self.param_grids = self._prepare_grids()

    # Preprocess and cache network topology for fast routing
    self._network_topology = self._preprocess_network_topology()

    # Pre-compute interpolation weights for meteorological forcing (station mode only)
    # Currently, only precipitation is used in the main loop
    if self.forcing_mode == "station":
        self._interpolation_weights = self._precompute_interpolation_weights(["precipitation"])
    # else: self._interpolation_weights already set to None above

    # Time indices cache (will be populated in run() when simulation period is known)
    self._time_indices_cache = None

    # Initialize state
    self.state = None

    logger.info(
        f"Simulation initialized: grid={self.nrows}x{self.ncols}, "
        f"dt={self.dt}s, network={len(self.network)} reaches"
    )

run(start_date, end_date, save_states_interval=None, save_report_interval=None)

Run simulation.

Parameters:

Name	Type	Description	Default
start_date	str | datetime	Simulation start date (YYYY-MM-DD or datetime)	required
end_date	str | datetime	Simulation end date (YYYY-MM-DD or datetime)	required
save_states_interval	int | None	Interval for saving states [s]. If None, use config value.	None
save_report_interval	int | None	Interval for saving reports [s]. If None, use config value.	None

Returns:

Type	Description
SimulationResults	SimulationResults object containing time series and final state

Source code in mobidic/core/simulation.py

def run(
    self,
    start_date: str | datetime,
    end_date: str | datetime,
    save_states_interval: int | None = None,
    save_report_interval: int | None = None,
) -> SimulationResults:
    """Run simulation.

    Args:
        start_date: Simulation start date (YYYY-MM-DD or datetime)
        end_date: Simulation end date (YYYY-MM-DD or datetime)
        save_states_interval: Interval for saving states [s]. If None, use config value.
        save_report_interval: Interval for saving reports [s]. If None, use config value.

    Returns:
        SimulationResults object containing time series and final state
    """

    logger.info("")
    logger.info("=" * 80)
    logger.info(f"MOBIDICpy v{__version__} - SIMULATION")
    logger.info("=" * 80)
    logger.info(f"Basin: {self.config.basin.id}")
    logger.info(f"Parameter set: {self.config.basin.paramset_id}")
    logger.info("")

    # Convert dates to datetime
    if isinstance(start_date, str):
        start_date = datetime.fromisoformat(start_date)
    if isinstance(end_date, str):
        end_date = datetime.fromisoformat(end_date)

    # Initialize state (skip if already set via set_initial_state())
    if self.state is None:
        self.state = self._initial_state()
    else:
        logger.info("Using pre-set initial state (skipping default initialization)")

    # Initialize results container
    results = SimulationResults(self.config, simulation=self)
    discharge_ts = []
    lateral_inflow_ts = []
    time_ts = []

    # Calculate number of time steps (inclusive of end_date)
    n_steps = int((end_date - start_date).total_seconds() / self.dt) + 1
    logger.info(f"Number of time steps: {n_steps}")

    # Pre-compute time indices for all simulation timesteps (station mode only)
    # Currently, only precipitation is used in the main loop
    simulation_times = pd.date_range(start=start_date, periods=n_steps, freq=f"{self.dt}s")
    if self.forcing_mode == "station":
        self._time_indices_cache = self._precompute_time_indices(simulation_times, variables=["precipitation"])
    else:
        self._time_indices_cache = None

    # Validate output_list if using list-based state output
    state_settings = self.config.output_states_settings
    if state_settings.output_states == "list" and state_settings.output_list:
        logger.debug("Validating output_list against simulation time steps")
        missing_states = []
        for date_str in state_settings.output_list:
            try:
                target_time = datetime.fromisoformat(date_str.replace(" ", "T"))
                # Check if target_time exists in simulation_times
                if target_time not in simulation_times:
                    missing_states.append(date_str)
            except (ValueError, TypeError) as e:
                logger.warning(f"Invalid datetime format in output_list: '{date_str}' - {e}")
                missing_states.append(date_str)

        if missing_states:
            logger.warning(
                f"The following {len(missing_states)} state(s) from output_list "
                f"will NOT be saved (not found among simulation time steps):"
            )
            for missing_state in missing_states:
                logger.warning(f"  - {missing_state}")
            logger.warning(
                f"Simulation runs from {start_date.isoformat()} to {end_date.isoformat()} with timestep={self.dt}s"
            )

    # Get total soil capacity
    wtot0 = self.wc0 + self.wg0

    # Cell area for unit conversion [m²]
    cell_area = self.resolution[0] * self.resolution[1]

    # Create ko mask: contributing pixels only (matching MATLAB mobidic_sid.m:224)
    # ko = find(isfinite(zz) & (ch>0));  % contributing pixels
    ko = np.isfinite(self.dtm) & (self.hillslope_reach_map >= 0)
    ko = np.where(ko.ravel("F"))[0]
    logger.debug(f"Contributing pixels (ko): {len(ko)} of {self.nrows * self.ncols} total cells")

    # Initialize flow variables for hillslope routing feedback (matching MATLAB mobidic_sid.m:1621-1625)
    # flr_prev: surface runoff from previous timestep [m/s] - will be routed to create pir
    # fld_prev: lateral flow from previous timestep [m/s] - will be routed to create pid
    flr_prev = np.zeros((self.nrows, self.ncols))
    fld_prev = np.zeros((self.nrows, self.ncols))

    # Calculate progress logging interval: either 20 steps total or every 30 seconds
    # Use whichever results in fewer log messages
    steps_per_intervals = max(1, n_steps // 20)
    progress_step_interval = steps_per_intervals
    progress_time_interval = 30.0  # seconds

    # Prepare states directory and initialize state writer (if enabled)
    state_settings = self.config.output_states_settings
    state_writer = None

    # Only instantiate StateWriter if state output is enabled
    if state_settings.output_states not in [None, "None"]:
        states_dir = Path(self.config.paths.states)
        states_dir.mkdir(parents=True, exist_ok=True)

        reservoir_size = len(self.reservoirs) if self.reservoirs is not None else 0
        state_writer = StateWriter(
            output_path=states_dir / "states.nc",
            grid_metadata=self.gisdata.metadata,
            network_size=len(self.network),
            output_states=self.config.output_states,
            flushing=state_settings.flushing,
            max_file_size=state_settings.max_file_size,
            add_metadata={
                "basin_id": self.config.basin.id,
                "paramset_id": self.config.basin.paramset_id,
            },
            reservoir_size=reservoir_size,
        )
        logger.info(f"State output enabled: {state_settings.output_states}")
    else:
        logger.info("State output disabled (output_states=None)")

    # Initialize meteo writer (if enabled)
    meteo_writer = None
    if self.config.output_forcing_data.meteo_data:
        output_dir = Path(self.config.paths.output)
        output_dir.mkdir(parents=True, exist_ok=True)

        from mobidic.io import MeteoWriter

        # Define which variables to save
        meteo_variables = ["precipitation"]

        meteo_writer = MeteoWriter(
            output_path=output_dir / "meteo_forcing.nc",
            grid_metadata=self.gisdata.metadata,
            variables=meteo_variables,
            add_metadata={
                "basin_id": self.config.basin.id,
                "paramset_id": self.config.basin.paramset_id,
                "start_date": start_date.isoformat(),
                "end_date": end_date.isoformat(),
            },
        )
        logger.info(f"Meteo data output enabled: {meteo_variables}")
    else:
        logger.info("Interpolated meteo data output disabled")

    # Main time loop
    logger.info("Starting simulation main loop")
    current_time = start_date
    simulation_start_time = time.time()
    last_log_time = simulation_start_time
    for step in range(n_steps):
        logger.debug(f"Time step {step + 1}/{n_steps}: {current_time}")

        # 1. Get meteorological forcing (interpolation for stations, direct sampling for rasters)
        try:
            precip = self._get_forcing_fn(
                current_time,
                "precipitation",
                weights_cache=self._interpolation_weights,
                time_step_index=step,
            )
        except (KeyError, ValueError):
            logger.warning(f"Precipitation data not found for {current_time}, using zero")
            precip = np.zeros((self.nrows, self.ncols))

        # Log precipitation statistics (precip is in m/s, convert to mm/h for readability)
        precip_mm_h = precip * 1000.0 * 3600.0  # m/s -> mm/h
        precip_valid = precip_mm_h[np.isfinite(precip_mm_h)]
        if len(precip_valid) > 0:
            logger.debug(
                f"Precipitation: mean={np.mean(precip_valid):.4f} mm/h, max={np.max(precip_valid):.4f} mm/h"
            )

        # 2. Calculate PET
        pet = self._calculate_pet(current_time)

        # 3. Save interpolated meteorological data (if enabled)
        if meteo_writer is not None:
            meteo_writer.append(current_time, precipitation=precip)

        # 4. Hillslope routing of previous timestep's flows (matching MATLAB mobidic_sid.m:1621-1625)
        # This must happen BEFORE soil mass balance to provide upstream contributions
        logger.debug("Routing previous timestep's flows through hillslope")

        # Convert to discharge for routing
        flr_prev_discharge = flr_prev * cell_area
        fld_prev_discharge = fld_prev * cell_area

        # Route through hillslope
        pir_discharge = hillslope_routing(flr_prev_discharge, self.flow_dir)
        pid_discharge = hillslope_routing(fld_prev_discharge, self.flow_dir)
        logger.debug("Hillslope routing completed")

        # Convert back to rates
        pir = pir_discharge / cell_area
        pid = pid_discharge / cell_area

        # 5. Soil water balance (cell-by-cell)
        logger.debug("Computing soil water balance")

        # Convert precipitation to depth over time step [m]
        precip_depth = precip * self.dt

        # Flatten 2D arrays to 1D (MATLAB: uses linear indexing with ko)
        # Extract only contributing pixels (ko) for processing (matching MATLAB mobidic_sid.m:733-736)
        wc_flat_full = self.state.wc.ravel("F")
        wg_flat_full = self.state.wg.ravel("F")
        wp_flat_full = self.state.wp.ravel("F") if self.state.wp is not None else None
        ws_flat_full = self.state.ws.ravel("F")

        # Extract ko cells only (matching MATLAB: Wc_ko = Wc(ko))
        wc_flat = wc_flat_full[ko]
        wg_flat = wg_flat_full[ko]
        wp_flat = wp_flat_full[ko] if wp_flat_full is not None else None
        ws_flat = ws_flat_full[ko]
        wc0_flat = self.wc0.ravel("F")[ko]
        wg0_flat = self.wg0.ravel("F")[ko]
        wtot0_flat = wtot0.ravel("F")[ko]
        precip_flat = precip_depth.ravel("F")[ko]
        pet_flat = pet.ravel("F")[ko] * self.dt
        # Multiply rate parameters by dt (matching MATLAB mobidic_sid.m:1681-1682)
        ks_flat = self.param_grids["ks"].ravel("F")[ko] * self.dt
        gamma_flat = self.param_grids["gamma"].ravel("F")[ko] * self.dt
        kappa_flat = self.param_grids["kappa"].ravel("F")[ko] * self.dt
        beta_flat = self.param_grids["beta"].ravel("F")[ko] * self.dt
        cha_flat = self.param_grids["cha"].ravel("F")[ko]
        f0_flat = self.param_grids["f0"].ravel("F")[ko]
        alpsur_flat = self.param_grids["alpsur"].ravel("F")[ko] * self.dt

        # Prepare routed flows from previous timestep (matching MATLAB mobidic_sid.m:1680)
        # Convert from [m/s] to [m] by multiplying by dt
        pir_flat = pir.ravel("F")[ko] * self.dt
        pid_flat = pid.ravel("F")[ko] * self.dt

        # Call soil_mass_balance with flattened arrays
        # Parameters have been pre-multiplied by dt above (lines 609-612, 616)
        (
            wc_out_flat,
            wg_out_flat,
            wp_out_flat,
            ws_out_flat,
            surface_runoff_flat,
            lateral_flow_depth_flat,
            et_flat,
            percolation_flat,
            capillary_flux_flat,
            wg_before_flat,
        ) = soil_mass_balance(
            wc=wc_flat,
            wc0=wc0_flat,
            wg=wg_flat,
            wg0=wg0_flat,
            wp=wp_flat,
            wp0=None,  # Plant reservoir disabled for now
            ws=ws_flat,
            ws0=None,
            wtot0=wtot0_flat,
            precipitation=precip_flat,
            surface_runoff_in=pir_flat,  # Routed surface runoff from previous timestep
            lateral_flow_in=pid_flat,  # Routed lateral flow from previous timestep
            potential_et=pet_flat,
            hydraulic_conductivity=ks_flat,
            hydraulic_conductivity_min=None,
            hydraulic_conductivity_max=None,
            channelized_fraction=cha_flat,
            surface_flow_exp=np.zeros_like(wc_flat),  # Not used
            lateral_flow_coeff=beta_flat,
            percolation_coeff=gamma_flat,
            absorption_coeff=kappa_flat,
            rainfall_fraction=f0_flat,
            et_shape=3.0,
            capillary_rise_enabled=False,
            test_mode=False,
            alpsur=alpsur_flat,
        )

        # Write ko results back to full grids (matching MATLAB: Wc(ko) = Wc_ko)
        wc_flat_full[ko] = wc_out_flat
        wg_flat_full[ko] = wg_out_flat
        if wp_out_flat is not None:
            wp_flat_full[ko] = wp_out_flat
        ws_flat_full[ko] = ws_out_flat

        # Initialize full output arrays for fluxes (preserve NaN for cells outside domain)
        surface_runoff_full = np.full(self.nrows * self.ncols, np.nan)
        lateral_flow_depth_full = np.full(self.nrows * self.ncols, np.nan)
        # Set non-ko cells within basin to zero
        surface_runoff_full[np.isfinite(self.dtm.ravel("F"))] = 0.0
        lateral_flow_depth_full[np.isfinite(self.dtm.ravel("F"))] = 0.0
        # Update ko cells with computed values
        surface_runoff_full[ko] = surface_runoff_flat
        lateral_flow_depth_full[ko] = lateral_flow_depth_flat

        # Reshape outputs back to 2D (Fortran order to match MATLAB)
        self.state.wc = wc_flat_full.reshape((self.nrows, self.ncols), order="F")
        self.state.wg = wg_flat_full.reshape((self.nrows, self.ncols), order="F")
        self.state.wp = (
            wp_flat_full.reshape((self.nrows, self.ncols), order="F") if wp_flat_full is not None else None
        )
        self.state.ws = ws_flat_full.reshape((self.nrows, self.ncols), order="F")
        surface_runoff = surface_runoff_full.reshape((self.nrows, self.ncols), order="F")
        lateral_flow_depth = lateral_flow_depth_full.reshape((self.nrows, self.ncols), order="F")

        # 5. Convert outputs from depth [m] to rate [m/s] (matching MATLAB mobidic_sid.m:1684)
        # flr: surface runoff rate [m/s] - analogous to MATLAB flr after division by dt
        # fld: lateral flow rate [m/s] - analogous to MATLAB fld after division by dt
        flr = surface_runoff / self.dt
        fld = lateral_flow_depth / self.dt

        # 6. Accumulate lateral inflow to reaches from surface runoff (matching MATLAB glob_route_day.m)
        # Convert to discharge for accumulation
        flr_discharge = flr * cell_area
        lateral_inflow = self._accumulate_lateral_inflow(flr_discharge)
        logger.debug("Lateral inflow accumulation completed")

        # Store lateral inflow in state
        self.state.lateral_inflow = lateral_inflow

        # Zero out flr for ALL cells that contributed to reaches (matching MATLAB glob_route_day.m line 33)
        # This prevents double-counting - flows from all contributing cells are consumed after accumulation
        flr[self.hillslope_reach_map >= 0] = 0.0
        # Note: fld is NOT zeroed - lateral flow continues to route between hillslope cells

        # Store flr and fld for next timestep's routing (at step 3)
        flr_prev = flr.copy()
        fld_prev = fld.copy()

        # 6. Reservoir routing (if reservoirs exist)
        # Prepare network topology for routing (may be modified if reservoirs exist)
        routing_network = self._network_topology

        if self.reservoirs is not None and self.state.reservoir_states is not None:
            logger.debug("Computing reservoir routing")

            # Call reservoir routing
            (
                self.state.reservoir_states,
                self.state.discharge,
                flr,
                fld,
                self.state.wg,
            ) = reservoir_routing(
                reservoirs_data=self.reservoirs.reservoirs,
                reservoir_states=self.state.reservoir_states,
                reach_discharge=self.state.discharge,
                surface_runoff=flr,
                lateral_flow=fld,
                soil_wg=self.state.wg,
                soil_wg0=self.wg0,
                current_time=current_time,
                dt=self.dt,
                cell_area=cell_area,
            )

            # Add reservoir outflows to lateral inflow of outlet reaches (matching MATLAB glob_route_day.m:46-53)
            for i, reservoir in enumerate(self.reservoirs.reservoirs):
                outlet_reach = reservoir.outlet_reach
                if outlet_reach is not None:
                    lateral_inflow[outlet_reach] += self.state.reservoir_states[i].outflow

            # Clear upstream connections for outlet reaches (matching MATLAB glob_route_day.m:48-50)
            # This prevents double-counting: the reservoir has already collected upstream volumes
            # Create modified topology with cleared upstream connections for outlet reaches
            routing_network = self._network_topology.copy()
            routing_network["upstream_1_idx"] = self._network_topology["upstream_1_idx"].copy()
            routing_network["upstream_2_idx"] = self._network_topology["upstream_2_idx"].copy()
            routing_network["n_upstream"] = self._network_topology["n_upstream"].copy()

            for reservoir in self.reservoirs.reservoirs:
                outlet_reach = reservoir.outlet_reach
                if outlet_reach is not None:
                    # Clear upstream connections (MATLAB: ret(j).ramimonte=[nan nan])
                    routing_network["upstream_1_idx"][outlet_reach] = -1
                    routing_network["upstream_2_idx"][outlet_reach] = -1
                    # Set upstream count to 0 (MATLAB: ram_fin(j) = 0)
                    routing_network["n_upstream"][outlet_reach] = 0

        # 7. Channel routing
        logger.debug("Computing channel routing")
        self.state.discharge, routing_state = linear_channel_routing(
            network=routing_network,
            discharge_initial=self.state.discharge,
            lateral_inflow=lateral_inflow,
            dt=self.dt,
        )

        # 8. Store results
        discharge_ts.append(self.state.discharge.copy())
        lateral_inflow_ts.append(lateral_inflow.copy())
        time_ts.append(current_time)

        # 9. Save intermediate states if configured
        if state_writer is not None and self._should_save_state(step, current_time):
            logger.debug(f"Appending state to buffer at step {step + 1}/{n_steps}")
            state_writer.append_state(self.state, current_time)

        # Log progress based on step interval or time interval (whichever comes first)
        current_wall_time = time.time()
        time_since_last_log = current_wall_time - last_log_time
        is_step_interval = (step + 1) % progress_step_interval == 0
        is_time_interval = time_since_last_log >= progress_time_interval
        is_final_step = step == n_steps - 1

        if is_step_interval or is_time_interval or is_final_step:
            progress_bar = _create_progress_bar(step + 1, n_steps, bar_length=20)
            q_mean = np.mean(self.state.discharge)
            q_max = np.max(self.state.discharge)
            date_str = current_time.strftime("%Y-%m-%d %H:%M")

            logger.info(
                f"{progress_bar} {step + 1}/{n_steps} | Simulation time: {date_str} | "
                f"Q_mean={q_mean:.3f} m³/s | Q_max={q_max:.3f} m³/s"
            )
            last_log_time = current_wall_time

        # Advance time
        current_time += timedelta(seconds=self.dt)

    # Store results
    results.time_series["discharge"] = np.array(discharge_ts)
    results.time_series["lateral_inflow"] = np.array(lateral_inflow_ts)
    results.time_series["time"] = time_ts
    results.final_state = self.state

    # Calculate elapsed time
    simulation_end_time = time.time()
    elapsed_seconds = simulation_end_time - simulation_start_time
    elapsed_minutes = elapsed_seconds / 60
    elapsed_hours = elapsed_minutes / 60

    # Format elapsed time message
    if elapsed_hours >= 1:
        elapsed_str = f"{elapsed_hours:.2f} hours"
    elif elapsed_minutes >= 1:
        elapsed_str = f"{elapsed_minutes:.2f} minutes"
    else:
        elapsed_str = f"{elapsed_seconds:.2f} seconds"

    logger.success(f"Simulation completed: {n_steps} time steps in {elapsed_str}")

    # Auto-save reports based on configuration
    output_dir = Path(self.config.paths.output)
    output_dir.mkdir(parents=True, exist_ok=True)

    # Determine reach selection
    settings = self.config.output_report_settings
    reach_selection = settings.reach_selection
    selected_reaches = None
    reach_file = None

    if settings.reach_selection == "file":
        reach_file = settings.sel_file
    elif settings.reach_selection == "list":
        selected_reaches = settings.sel_list

    # Get output format from configuration
    output_format = self.config.output_report_settings.output_format
    file_extension = "csv" if output_format.lower() == "csv" else "parquet"

    # Save discharge report if enabled
    if self.config.output_report.discharge:
        start_str = start_date.strftime("%Y%m%d")
        end_str = end_date.strftime("%Y%m%d")
        discharge_path = output_dir / f"discharge_{start_str}_{end_str}.{file_extension}"
        logger.info("Exporting discharges")
        results.save_report(
            output_path=discharge_path,
            reach_selection=reach_selection,
            selected_reaches=selected_reaches,
            reach_file=reach_file,
            output_format=output_format,
        )

    # Save lateral inflow report if enabled
    if self.config.output_report.lateral_inflow:
        lateral_inflow_path = output_dir / f"lateral_inflow_{start_str}_{end_str}.{file_extension}"
        logger.info("Exporting lateral inflows")
        results.save_lateral_inflow_report(
            output_path=lateral_inflow_path,
            reach_selection=reach_selection,
            selected_reaches=selected_reaches,
            reach_file=reach_file,
            output_format=output_format,
        )

    # Save final state if enabled (for "final" mode, only save the last state)
    if state_writer is not None:
        if state_settings.output_states == "final":
            final_time = results.time_series["time"][-1]
            logger.info("Saving final state")
            state_writer.append_state(self.state, final_time)

        # Close the state writer (flushes any remaining buffered states)
        state_writer.close()

    # Close the meteo writer (writes all buffered data to NetCDF)
    if meteo_writer is not None:
        meteo_writer.close()

    return results

set_initial_state(state=None, state_file=None, time_index=-1)

Set the initial simulation state from a previous simulation.

This method allows restarting a simulation from a previously saved state, enabling warm starts, multi-stage simulations, or continuation after interruption.

Parameters:

Name	Type	Description	Default
state	SimulationState | None	SimulationState object to use as initial state. If provided, state_file is ignored.	None
state_file	str | Path | None	Path to NetCDF state file to load. Used if state is None.	None
time_index	int	Time index to load from state file (default: -1 for last timestep). Only used when loading from state_file.	-1

Raises:

Type	Description
ValueError	If neither state nor state_file is provided, or if state_file doesn’t exist.

Examples:

>>> # Method 1: Set state directly from SimulationState object
>>> sim.set_initial_state(state=previous_state)
>>>
>>> # Method 2: Load from state file (last timestep)
>>> sim.set_initial_state(state_file="states.nc")
>>>
>>> # Method 3: Load specific timestep from state file
>>> sim.set_initial_state(state_file="states.nc", time_index=10)

Source code in mobidic/core/simulation.py

def set_initial_state(
    self,
    state: SimulationState | None = None,
    state_file: str | Path | None = None,
    time_index: int = -1,
) -> None:
    """Set the initial simulation state from a previous simulation.

    This method allows restarting a simulation from a previously saved state,
    enabling warm starts, multi-stage simulations, or continuation after interruption.

    Args:
        state: SimulationState object to use as initial state. If provided, state_file is ignored.
        state_file: Path to NetCDF state file to load. Used if state is None.
        time_index: Time index to load from state file (default: -1 for last timestep).
            Only used when loading from state_file.

    Raises:
        ValueError: If neither state nor state_file is provided, or if state_file doesn't exist.

    Examples:
        >>> # Method 1: Set state directly from SimulationState object
        >>> sim.set_initial_state(state=previous_state)
        >>>
        >>> # Method 2: Load from state file (last timestep)
        >>> sim.set_initial_state(state_file="states.nc")
        >>>
        >>> # Method 3: Load specific timestep from state file
        >>> sim.set_initial_state(state_file="states.nc", time_index=10)
    """
    if state is not None:
        # Use provided state directly
        logger.info("Setting initial state from provided SimulationState object")
        self.state = state
        logger.success(
            f"Initial state set. State contains: "
            f"Wc={np.nanmean(state.wc) * 1000:.1f} mm, "
            f"Wg={np.nanmean(state.wg) * 1000:.1f} mm, "
            f"Ws={np.nanmean(state.ws) * 1000:.1f} mm, "
            f"Q_mean={np.mean(state.discharge):.3f} m³/s"
        )
    elif state_file is not None:
        # Load from state file
        state_path = Path(state_file)
        if not state_path.exists():
            raise ValueError(f"State file not found: {state_path}")

        logger.info(f"Loading initial state from file: {state_path}")

        from mobidic.io import load_state

        loaded_state, state_time, metadata = load_state(
            input_path=state_path,
            network_size=len(self.network),
            time_index=time_index,
            config=self.config,
            gisdata=self.gisdata,
        )

        self.state = loaded_state
        logger.success(
            f"Initial state loaded from {state_path} at time {state_time}. "
            f"State contains: "
            f"Wc={np.nanmean(loaded_state.wc) * 1000:.1f} mm, "
            f"Wg={np.nanmean(loaded_state.wg) * 1000:.1f} mm, "
            f"Ws={np.nanmean(loaded_state.ws) * 1000:.1f} mm, "
            f"Q_mean={np.mean(loaded_state.discharge):.3f} m³/s"
        )
    else:
        raise ValueError("Either 'state' or 'state_file' must be provided")

Container for simulation state variables.

Source code in mobidic/core/simulation.py

class SimulationState:
    """Container for simulation state variables."""

    def __init__(
        self,
        wc: np.ndarray,
        wg: np.ndarray,
        wp: np.ndarray | None,
        ws: np.ndarray,
        discharge: np.ndarray,
        lateral_inflow: np.ndarray,
        reservoir_states: list | None = None,
    ):
        """Initialize simulation state.

        Args:
            wc: Capillary water content [m]
            wg: Gravitational water content [m]
            wp: Plant/canopy water content [m] (None to disable)
            ws: Surface water content [m]
            discharge: River discharge for each reach [m³/s]
            lateral_inflow: Lateral inflow to each reach [m³/s]
            reservoir_states: List of ReservoirState objects (None if no reservoirs)
        """
        self.wc = wc
        self.wg = wg
        self.wp = wp
        self.ws = ws
        self.discharge = discharge
        self.lateral_inflow = lateral_inflow
        self.reservoir_states = reservoir_states

__init__(wc, wg, wp, ws, discharge, lateral_inflow, reservoir_states=None)

Initialize simulation state.

Parameters:

Name	Type	Description	Default
wc	ndarray	Capillary water content [m]	required
wg	ndarray	Gravitational water content [m]	required
wp	ndarray | None	Plant/canopy water content [m] (None to disable)	required
ws	ndarray	Surface water content [m]	required
discharge	ndarray	River discharge for each reach [m³/s]	required
lateral_inflow	ndarray	Lateral inflow to each reach [m³/s]	required
reservoir_states	list | None	List of ReservoirState objects (None if no reservoirs)	None

Source code in mobidic/core/simulation.py

def __init__(
    self,
    wc: np.ndarray,
    wg: np.ndarray,
    wp: np.ndarray | None,
    ws: np.ndarray,
    discharge: np.ndarray,
    lateral_inflow: np.ndarray,
    reservoir_states: list | None = None,
):
    """Initialize simulation state.

    Args:
        wc: Capillary water content [m]
        wg: Gravitational water content [m]
        wp: Plant/canopy water content [m] (None to disable)
        ws: Surface water content [m]
        discharge: River discharge for each reach [m³/s]
        lateral_inflow: Lateral inflow to each reach [m³/s]
        reservoir_states: List of ReservoirState objects (None if no reservoirs)
    """
    self.wc = wc
    self.wg = wg
    self.wp = wp
    self.ws = ws
    self.discharge = discharge
    self.lateral_inflow = lateral_inflow
    self.reservoir_states = reservoir_states

Container for simulation results.

Source code in mobidic/core/simulation.py

class SimulationResults:
    """Container for simulation results."""

    def __init__(self, config: MOBIDICConfig, simulation: "Simulation | None" = None):
        """Initialize results container.

        Args:
            config: MOBIDIC configuration
            simulation: Simulation object (needed for saving states/reports)
        """
        self.config = config
        self.simulation = simulation
        self.time_series = {}  # Store time series data
        self.final_state = None  # Store final state

    def save_report(
        self,
        output_path: str | Path,
        reach_selection: str = "all",
        selected_reaches: list[int] | None = None,
        reach_file: str | Path | None = None,
        add_metadata: dict[str, Any] | None = None,
        output_format: str = "Parquet",
    ) -> None:
        """Save discharge time series to file (Parquet or CSV).

        Args:
            output_path: Path to output file
            reach_selection: "all", "file", or "list"
            selected_reaches: List of reach IDs (if reach_selection="list")
            reach_file: Path to JSON file containing reach IDs (if reach_selection="file")
            add_metadata: Additional metadata to save (optional)
            output_format: Output format: "Parquet" or "csv" (default: "Parquet")
        """
        if "discharge" not in self.time_series:
            raise ValueError("No discharge data to save. Run simulation first.")

        if self.simulation is None:
            raise ValueError("Cannot save report without simulation object")

        from mobidic.io import save_discharge_report

        save_discharge_report(
            discharge_timeseries=self.time_series["discharge"],
            time_stamps=self.time_series["time"],
            network=self.simulation.network,
            output_path=output_path,
            reach_selection=reach_selection,
            selected_reaches=selected_reaches,
            reach_file=reach_file,
            add_metadata=add_metadata,
            output_format=output_format,
        )

    def save_lateral_inflow_report(
        self,
        output_path: str | Path,
        reach_selection: str = "all",
        selected_reaches: list[int] | None = None,
        reach_file: str | Path | None = None,
        output_format: str = "Parquet",
    ) -> None:
        """Save lateral inflow time series to file (Parquet or CSV).

        Args:
            output_path: Path to output file
            reach_selection: "all", "file", or "list"
            selected_reaches: List of reach IDs (if reach_selection="list")
            reach_file: Path to JSON file containing reach IDs (if reach_selection="file")
            output_format: Output format: "Parquet" or "csv" (default: "Parquet")
        """
        if "lateral_inflow" not in self.time_series:
            raise ValueError("No lateral inflow data to save. Run simulation first.")

        if self.simulation is None:
            raise ValueError("Cannot save lateral inflow report without simulation object")

        from mobidic.io import save_lateral_inflow_report

        save_lateral_inflow_report(
            lateral_inflow_timeseries=self.time_series["lateral_inflow"],
            time_stamps=self.time_series["time"],
            network=self.simulation.network,
            output_path=output_path,
            reach_selection=reach_selection,
            selected_reaches=selected_reaches,
            reach_file=reach_file,
            output_format=output_format,
        )

__init__(config, simulation=None)

Initialize results container.

Parameters:

Name	Type	Description	Default
config	MOBIDICConfig	MOBIDIC configuration	required
simulation	Simulation | None	Simulation object (needed for saving states/reports)	None

Source code in mobidic/core/simulation.py

def __init__(self, config: MOBIDICConfig, simulation: "Simulation | None" = None):
    """Initialize results container.

    Args:
        config: MOBIDIC configuration
        simulation: Simulation object (needed for saving states/reports)
    """
    self.config = config
    self.simulation = simulation
    self.time_series = {}  # Store time series data
    self.final_state = None  # Store final state

save_lateral_inflow_report(output_path, reach_selection='all', selected_reaches=None, reach_file=None, output_format='Parquet')

Save lateral inflow time series to file (Parquet or CSV).

Parameters:

Name	Type	Description	Default
output_path	str | Path	Path to output file	required
reach_selection	str	“all”, “file”, or “list”	'all'
selected_reaches	list[int] | None	List of reach IDs (if reach_selection=”list”)	None
reach_file	str | Path | None	Path to JSON file containing reach IDs (if reach_selection=”file”)	None
output_format	str	Output format: “Parquet” or “csv” (default: “Parquet”)	'Parquet'

Source code in mobidic/core/simulation.py

def save_lateral_inflow_report(
    self,
    output_path: str | Path,
    reach_selection: str = "all",
    selected_reaches: list[int] | None = None,
    reach_file: str | Path | None = None,
    output_format: str = "Parquet",
) -> None:
    """Save lateral inflow time series to file (Parquet or CSV).

    Args:
        output_path: Path to output file
        reach_selection: "all", "file", or "list"
        selected_reaches: List of reach IDs (if reach_selection="list")
        reach_file: Path to JSON file containing reach IDs (if reach_selection="file")
        output_format: Output format: "Parquet" or "csv" (default: "Parquet")
    """
    if "lateral_inflow" not in self.time_series:
        raise ValueError("No lateral inflow data to save. Run simulation first.")

    if self.simulation is None:
        raise ValueError("Cannot save lateral inflow report without simulation object")

    from mobidic.io import save_lateral_inflow_report

    save_lateral_inflow_report(
        lateral_inflow_timeseries=self.time_series["lateral_inflow"],
        time_stamps=self.time_series["time"],
        network=self.simulation.network,
        output_path=output_path,
        reach_selection=reach_selection,
        selected_reaches=selected_reaches,
        reach_file=reach_file,
        output_format=output_format,
    )

save_report(output_path, reach_selection='all', selected_reaches=None, reach_file=None, add_metadata=None, output_format='Parquet')

Save discharge time series to file (Parquet or CSV).

Parameters:

Name	Type	Description	Default
output_path	str | Path	Path to output file	required
reach_selection	str	“all”, “file”, or “list”	'all'
selected_reaches	list[int] | None	List of reach IDs (if reach_selection=”list”)	None
reach_file	str | Path | None	Path to JSON file containing reach IDs (if reach_selection=”file”)	None
add_metadata	dict[str, Any] | None	Additional metadata to save (optional)	None
output_format	str	Output format: “Parquet” or “csv” (default: “Parquet”)	'Parquet'

Source code in mobidic/core/simulation.py

def save_report(
    self,
    output_path: str | Path,
    reach_selection: str = "all",
    selected_reaches: list[int] | None = None,
    reach_file: str | Path | None = None,
    add_metadata: dict[str, Any] | None = None,
    output_format: str = "Parquet",
) -> None:
    """Save discharge time series to file (Parquet or CSV).

    Args:
        output_path: Path to output file
        reach_selection: "all", "file", or "list"
        selected_reaches: List of reach IDs (if reach_selection="list")
        reach_file: Path to JSON file containing reach IDs (if reach_selection="file")
        add_metadata: Additional metadata to save (optional)
        output_format: Output format: "Parquet" or "csv" (default: "Parquet")
    """
    if "discharge" not in self.time_series:
        raise ValueError("No discharge data to save. Run simulation first.")

    if self.simulation is None:
        raise ValueError("Cannot save report without simulation object")

    from mobidic.io import save_discharge_report

    save_discharge_report(
        discharge_timeseries=self.time_series["discharge"],
        time_stamps=self.time_series["time"],
        network=self.simulation.network,
        output_path=output_path,
        reach_selection=reach_selection,
        selected_reaches=selected_reaches,
        reach_file=reach_file,
        add_metadata=add_metadata,
        output_format=output_format,
    )

Simulation loop

The main simulation loop performs the following operations for each time step:

1. Get forcing: Precipitation from station data (interpolated to grid using IDW/nearest) or from raster data (direct sampling)
2. Calculate PET: Simple 1 mm/day method (constant rate)
3. Save interpolated meteo (optional): Export interpolated grids when using station-based forcing
4. Route previous flows: Hillslope routing of surface runoff and lateral flow from previous timestep
5. Soil water balance: Four-reservoir hillslope water balance with routed inflows
6. Reservoir routing (if configured): Update reservoir volumes, calculate regulated discharge, zero basin fluxes
7. Accumulate to reaches: Accumulate surface runoff contributions to river reaches
8. Channel routing: Linear reservoir routing through river network
9. Store results: Save discharge and lateral inflow time series
10. Output states: Optionally save states (with automatic chunking if needed)
11. Update and advance: Store flow fields for next timestep and advance simulation time

Key feature: The simulation uses a feedback loop where flows from timestep t are routed through the hillslope at timestep t+1 before entering the soil water balance. This ensures proper spatial connectivity of overland flow.

Performance

• Meteorological interpolation caching: Pre-computes time indices and spatial weights for all timesteps
• Numba JIT compilation: Hillslope and channel routing use compiled kernels
• Memory efficiency: State variables use NumPy arrays with F-contiguous memory layout
• Contributing pixels optimization: Processes only cells that can contribute flow to river network
• Network topology caching: Pre-extracts network structure to numpy arrays for fast routing
• Progress logging: Adaptive logging interval (max 20 logs or every 30s) with text-based progress bar
• Automatic file chunking: State files automatically split when reaching size limit

Examples

Basic simulation with station-based forcing

from mobidic import load_config, load_gisdata, Simulation, MeteoData

# Load configuration and data
config = load_config("config.yaml")
gisdata = load_gisdata("gisdata.nc", "network.parquet")
forcing = MeteoData.from_netcdf("meteo.nc")

# Create simulation
sim = Simulation(gisdata, forcing, config)

# Run simulation
results = sim.run("2020-01-01", "2020-12-31")

# Save discharge report
results.save_report("discharge.parquet")

# Save lateral inflow report
results.save_lateral_inflow_report("lateral_inflow.parquet")

Simulation with raster-based forcing

from mobidic import load_config, load_gisdata, Simulation, MeteoRaster

# Load configuration and data
config = load_config("config.yaml")
gisdata = load_gisdata("gisdata.nc", "network.parquet")

# Load raster forcing (preload into memory for fast access)
forcing = MeteoRaster.from_netcdf("meteo_raster.nc")

# Create simulation (automatically detects raster mode)
sim = Simulation(gisdata, forcing, config)

# Run simulation (no interpolation needed, uses direct sampling)
results = sim.run("2020-01-01", "2020-12-31")

# Save results
results.save_report("discharge.parquet")

Export and use interpolated meteorological data

from mobidic import MeteoData, MeteoRaster, Simulation

# Run 1: Station-based with forcing output enabled
config.output_forcing_data.meteo_data = True
forcing_stations = MeteoData.from_netcdf("meteo_stations.nc")
sim1 = Simulation(gisdata, forcing_stations, config)
results1 = sim1.run("2020-01-01", "2020-12-31")
# Forcing data saved to: output/meteo_forcing.nc

# Run 2: Use exported raster forcing (faster, identical results)
config.output_forcing_data.meteo_data = False
forcing_raster = MeteoRaster.from_netcdf("output/meteo_forcing.nc")
sim2 = Simulation(gisdata, forcing_raster, config)
results2 = sim2.run("2020-01-01", "2020-12-31")

Warm start (resume from saved simulation state)

The simulation supports warm start capability, allowing you to resume from previously saved states. This is useful for:

• Multi-stage simulations: Spin-up period followed by analysis period
• Interrupted simulations: Resume after crashes or timeouts
• Ensemble runs: Start multiple simulations from calibrated initial states
• Seasonal forecasts: Initialize from observed states

from mobidic import Simulation, load_state

# Method 1: Load state from file and set before running
sim = Simulation(gisdata, forcing, config)
sim.set_initial_state(state_file="spinup_states.nc", time_index=-1)  # Use last timestep
results = sim.run("2020-06-01", "2020-12-31")

# Method 2: Load state object and use directly
from mobidic.io import load_state
state, time, metadata = load_state("spinup_states.nc", network_size=1235)
sim.set_initial_state(state=state)
results = sim.run("2020-06-01", "2020-12-31")

# Method 3: Multi-stage simulation
# Stage 1: Spin-up (1 year)
sim1 = Simulation(gisdata, forcing, config)
results1 = sim1.run("2019-01-01", "2019-12-31")  # Saves states.nc

# Stage 2: Analysis period (resume from spin-up)
sim2 = Simulation(gisdata, forcing, config)
sim2.set_initial_state(state_file="output/states.nc")  # Load last state from spin-up
results2 = sim2.run("2020-01-01", "2020-12-31")

Working with large simulations and chunked states

For long simulations that generate a large number of states, the simulation automatically creates chunked files:

from mobidic import Simulation

# Configure for large simulation
# Edit config.yaml:
# output_states_settings:
#   output_states: "all"
#   flushing: 100          # Flush every 100 timesteps (required for chunking)
#   max_file_size: 500.0   # Create new chunk at 500 MB
#   output_interval: 3600  # Save state every hour

config = load_config("config.yaml")
sim = Simulation(gisdata, forcing, config)

# Run long simulation (e.g., 10 years at 15-minute timesteps)
results = sim.run("2010-01-01", "2020-12-31")

# This creates multiple chunk files:
# - output/states_001.nc (500 MB)
# - output/states_002.nc (500 MB)
# - output/states_003.nc (350 MB)

# Resume from last chunk (automatically detected)
sim2 = Simulation(gisdata, forcing, config)
sim2.set_initial_state(state_file="output/states.nc")  # Auto-finds states_001.nc
results2 = sim2.run("2021-01-01", "2021-12-31")

Custom report selection

Control which reaches to include in output reports:

# Method 1: Save all reaches
results.save_report(
    "discharge_all.parquet",
    reach_selection="all"
)

# Method 2: Save reaches from file
results.save_report(
    "discharge_selected.parquet",
    reach_selection="file",
    reach_file="reach_ids.json"
)

# Method 3: Save specific reaches by mobidic_id
selected_reaches = [0, 10, 25, 100, 500]  # mobidic_id values
results.save_report(
    "discharge_selected.parquet",
    reach_selection="list",
    selected_reaches=selected_reaches
)

# Method 4: Export as CSV instead of Parquet
results.save_report(
    "discharge.csv",
    reach_selection="all",
    output_format="csv"
)

Configuration

The simulation behavior is controlled by the configuration file. Key sections:

Output states configuration

output_states:
  Wc: true              # Save capillary water
  Wg: true              # Save gravitational water
  Wp: false             # Plant water (not yet implemented)
  Ws: true              # Save surface water
  discharge: true       # Save channel discharge
  lateral_inflow: true  # Save lateral inflow to reaches
  reservoir_states: true  # Save reservoir states (if reservoirs configured)

output_states_settings:
  output_format: "netCDF"
  output_states: "final"  # Options: "final", "all", "list", "None"
  flushing: 10            # Flush to disk every N timesteps (-1 = only at end)
  max_file_size: 500.0    # Maximum file size in MB (chunking threshold)
  output_interval: 3600   # Save interval in seconds (for "all" mode)
  output_list:            # List of specific datetimes (for "list" mode)
    - "2020-06-01 00:00:00"
    - "2020-12-31 23:45:00"

Important notes about state output: - Chunking requires flushing > 0: Set flushing to a positive value (e.g., 10, 50, 100) for file chunking to work effectively - flushing=-1 disables chunking: All data is written at the end in one operation, preventing chunking - File size may slightly exceed limit: Files can exceed max_file_size by up to one flush worth of data

Output reports configuration

output_report:
  discharge: true        # Export discharge time series
  lateral_inflow: true   # Export lateral inflow time series

output_report_settings:
  output_format: "Parquet"     # Options: "Parquet", "csv"
  reach_selection: "all"       # Options: "all", "file", "list"
  sel_list: [0, 10, 25, 100]  # List of reach IDs (for "list" mode)
  sel_file: "reaches.json"     # Path to file with reach IDs (for "file" mode)

Reservoir configuration (optional)

parameters:
  reservoirs:
    res_shape: reservoirs/reservoirs.shp
    stage_storage: reservoirs/stage_storage.csv
    regulation_curves: reservoirs/regulation_curves.csv
    regulation_schedule: reservoirs/regulation_schedule.csv

initial_conditions:
  reservoir_volumes: reservoirs/initial_volumes.csv  # Optional

paths:
  reservoirs: output/reservoirs.parquet  # Consolidated reservoir data

When reservoirs are configured: - Reservoir polygons are rasterized to identify basin pixels - Surface runoff and lateral flow are zeroed in reservoir basins - Total inflow is computed from upstream discharge and basin contributions - Reservoir volume is updated based on mass balance - Stage is calculated from volume using cubic spline interpolation - Regulated discharge is determined from time-varying stage-discharge curves - Reservoir outflow is added to outlet reach lateral inflow - Inlet reach discharge is zeroed to prevent double-counting

Simulation configuration

simulation:
  timestep: 900                      # Time step in seconds (15 minutes)
  precipitation_interp: "IDW"        # Options: "IDW", "Nearest" (for station-based forcing)

output_forcing_data:
  meteo_data: false                  # Save meteorological forcing grids
                                     # Output file: {output_dir}/meteo_forcing.nc

Notes: - precipitation_interp only applies when using station-based forcing (MeteoData) - When using raster-based forcing (MeteoRaster), interpolation is skipped entirely

Implemented modules

• Preprocessing - GIS data and reservoir preprocessing
• Soil Water Balance - Hillslope water balance
• Routing - Hillslope, channel, and reservoir routing
• State I/O - NetCDF state export/import
• Report I/O - Time series export/import