Energy balance

The energy balance module computes the surface energy budget at the land surface and calculates the potential evapotranspiration (PET) before the soil-water balance calculations. Surface (\(T_s\)) and deep-soil (\(T_d\)) temperatures are calculated as additional state variables.

Overview

The currently available scheme is the 1-layer (1L) analytical formulation. It solves a linearised Fourier decomposition of the surface energy budget, with the diurnal cycle of incoming shortwave radiation represented as a sinusoid plus a constant term. Each model timestep is split at the sunrise and sunset boundaries so that day and night sub-periods are integrated separately (radiation forcing is only present in the day sub-period).

The energy budget includes:

Net shortwave radiation (corrected by surface albedo \(\alpha\))
Net longwave radiation (Stefan-Boltzmann linearised around mean air temperature)
Sensible heat flux (turbulent exchange coefficient \(C_H\) times wind speed)
Latent heat flux (Magnus formula for saturation specific humidity)
Ground heat flux (1-D conduction with deep-temperature boundary \(T_{const}\))

Two-step calculation

For each timestep the energy balance is solved in two steps:

Initial (soil saturation assumption): calculates the energy balance assuming no water limitation. The ratio \(\eta = \text{ET} / \text{PET}\) is set to 1.
Second step (after the soil water balance): re-computes \(T_s\) and \(T_d\) using the actual ratio \(\eta\) obtained from the soil module.

The deep-temperature value evaluated at sunrise is preserved across timesteps and is used as the starting point for the day sub-period in the second step calculation.

Diurnal radiation cycle

Daily incoming shortwave radiation \(R_s\) is decomposed into:

\[ R_s(t) = A \cdot \sin(\omega t + \varphi) + C \]

with \(\omega = 2\pi / 86400\) rad/s. Two modes are supported:

average: \(R_s\) is the day-average value; the amplitude \(A\) and constant \(C\) are fitted so that the integral between sunrise and sunset returns the daily mean.
instant: \(R_s\) is interpreted as an instantaneous value at sunrise time.

Sunrise and sunset hours are obtained by bisection on the solar elevation, computed from the basin baricenter latitude/longitude and the Julian day.

Land-cover adjustment via Kc

When a Corine Land Cover raster (raster_files.CLC) is provided, the turbulent exchange coefficient \(C_H\) is scaled by the monthly FAO crop coefficient \(K_c\) before the energy balance solver is called:

\[C_H^{\text{adj}} = K_c \cdot C_H\]

\(C_H^{\text{adj}}\) is then used in the calculation of both sensible heat flux \(H\) and the latent heat flux \(LE\). See Crop coefficients (Kc) for the full description and configuration options.

Configuration

The energy balance is enabled under simulation.energy_balance and is configured through parameters.energy, the optional CH / Alb rasters, and the basin.baricenter coordinates required to compute solar hours.

basin:
  baricenter:
    lon: 11.10   # Longitude of basin baricenter [deg. East]
    lat: 43.77   # Latitude of basin baricenter [deg. North]

simulation:
  # Energy balance scheme [None | 1L]
  energy_balance: 1L

raster_files:
  # OPTIONAL: Grid of turbulent exchange coefficient for heat, non dimensional
  CH: example/raster/CH.tif

  # OPTIONAL: Grid of surface albedo, non dimensional
  Alb: example/raster/Alb.tif

parameters:
  energy:
    Tconst: 290.0     # Deep ground temperature [K] (default: 290.0)
    kaps:   2.5       # Soil thermal conductivity [W/m/K] (default: 2.5)
    nis:    0.8e-6    # Soil thermal diffusivity [m^2/s] (default: 0.8e-6)
    CH:     1.0e-3    # Default turbulent exchange coefficient (used if CH raster missing)
    Alb:    0.2       # Default surface albedo (used if Alb raster missing)

  multipliers:
    CH_factor: 1.0    # Calibration multiplier applied to CH raster/scalar

output_states:
  surface_temperature: true   # Save Ts grid to state file
  ground_temperature:  true   # Save Td grid to state file

output_forcing_data:
  meteo_data: true            # Saves precipitation, temperature_min/max, humidity,
                              # wind_speed, radiation, and pet_c (Kc-adjusted PET) to meteo_forcing.nc

Required meteorological forcing. When simulation.energy_balance == "1L", the simulation needs the following variables in addition to precipitation:

temperature_min, temperature_max — daily air temperature extremes [°C]
humidity — relative humidity [%]
wind_speed — wind speed [m/s]
radiation — incoming shortwave radiation [W/m²]

Both station-based (MeteoData) and raster-based (MeteoRaster) inputs are supported.

Integration with the simulation loop

When simulation.energy_balance == "1L":

The simulation initialises \(T_s = T_d = T_{air,lin}\) on the first step (when no warm start is used) or restores them from the loaded state file. td_rise is initialised from \(T_d\).
Each timestep, the initial step is run with \(\eta = 1\) to obtain PET and a tentative \((T_s, T_d)\).
PET values are fed into the soil water balance, which produces the actual evapotranspiration ET.
The second step is run with \(\eta = \text{ET} / \text{PET}\), restarting the day sub-period from the initial td_rise. The result overwrites \((T_s, T_d)\) only on water-limited cells.
\(T_s\) and \(T_d\) are written to the state file when the corresponding output_states flags are enabled.

ET/PET when using raster forcing, to speed up simulations

When the input MeteoRaster contains either pet or pet_c variable, the simulation skips the energy balance entirely regardless of the simulation.energy_balance setting, and reads the evapotranspiration values directly from the raster. In this case, the simulation speeds up significantly. In this mode, the temperature/humidity/wind/radiation variables are not required in the input file. See Crop coefficients (Kc) for details.

Model status

None — energy balance disabled, constant 1 mm/day PET is used
1L — single-layer analytical Fourier scheme (this module)
5L — multi-layer soil temperature profile (not yet implemented)
Snow — snow accumulation and melt (not yet implemented)

Functions

`mobidic.core.energy_balance.compute_energy_balance_1l(ts, td, td_rise, rs, u, tair_max, tair_min, qair, ch, alb, kaps, nis, tcost, pair, ctim_s, ftim_s, hrise_s, hset_s, etrsuetp, dt, reentry=False)`

Run the 1L surface energy balance for a single model timestep.

Day and night sub-periods are solved separately (day uses radiation forcing, night does not).

Parameters:

Name	Type	Description	Default
`ts`	`NDArray[float64]`	Surface temperature [K] at the beginning of the step (per cell).	required
`td`	`NDArray[float64]`	Deep-soil temperature [K] at the beginning of the step (per cell).	required
`td_rise`	`NDArray[float64]`	`td` state evaluated at sunrise [K]. Used as the starting deep temperature for the day-period re-computation. In the initial (`reentry=False`) pass the function updates this field; in re-entry it is read-only.	required
`rs`	`NDArray[float64]`	Incoming shortwave radiation [W/m²] (per cell).	required
`u`	`NDArray[float64]`	Wind speed [m/s] (per cell).	required
`tair_max`	`NDArray[float64]`	Maximum air temperature [K] (per cell).	required
`tair_min`	`NDArray[float64]`	Minimum air temperature [K] (per cell).	required
`qair`	`NDArray[float64]`	Relative humidity [0-1] (per cell).	required
`ch`	`NDArray[float64]`	Turbulent exchange coefficient [-] (per cell).	required
`alb`	`NDArray[float64]`	Surface albedo [-] (per cell).	required
`kaps`	`float`	Soil thermal conductivity [W/m/K].	required
`nis`	`float`	Soil thermal diffusivity [m²/s].	required
`tcost`	`float`	Deep-soil constant temperature [K].	required
`pair`	`float`	Air pressure [mb].	required
`ctim_s`	`float`	Current-time seconds of day [s] (start of step).	required
`ftim_s`	`float`	Future-time seconds of day [s] (end of step = ctim_s + dt).	required
`hrise_s`	`float`	Sunrise time [s of day].	required
`hset_s`	`float`	Sunset time [s of day].	required
`etrsuetp`	`float \| NDArray[float64]`	Water-to-energy ET ratio [-]. `1.0` for the preliminary saturated pass; actual `ETr/ETp` for re-entry.	required
`dt`	`float`	Model timestep [s].	required
`reentry`	`bool`	If True, skip the pre-sunrise night sub-period (physics is independent of `etrsuetp`) and start day computation from `td_rise` rather than the current `td`.	`False`

Returns:

Name	Type	Description
`ts_new`	`NDArray[float64]`	Updated surface temperature [K].
`td_new`	`NDArray[float64]`	Updated deep temperature [K].
`etp`	`NDArray[float64]`	Potential evapotranspiration over the timestep [m].
`td_rise_new`	`NDArray[float64]`	Updated `td` at sunrise (for later re-entry).

Source code in mobidic/core/energy_balance.py

def compute_energy_balance_1l(
    ts: NDArray[np.float64],
    td: NDArray[np.float64],
    td_rise: NDArray[np.float64],
    rs: NDArray[np.float64],
    u: NDArray[np.float64],
    tair_max: NDArray[np.float64],
    tair_min: NDArray[np.float64],
    qair: NDArray[np.float64],
    ch: NDArray[np.float64],
    alb: NDArray[np.float64],
    kaps: float,
    nis: float,
    tcost: float,
    pair: float,
    ctim_s: float,
    ftim_s: float,
    hrise_s: float,
    hset_s: float,
    etrsuetp: float | NDArray[np.float64],
    dt: float,
    reentry: bool = False,
) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
    """Run the 1L surface energy balance for a single model timestep.

    Day and night sub-periods are solved separately
    (day uses radiation forcing, night does not).

    Args:
        ts: Surface temperature [K] at the beginning of the step (per cell).
        td: Deep-soil temperature [K] at the beginning of the step (per cell).
        td_rise: ``td`` state evaluated at sunrise [K]. Used as the starting
            deep temperature for the day-period re-computation. In the
            initial (``reentry=False``) pass the function updates this field;
            in re-entry it is read-only.
        rs: Incoming shortwave radiation [W/m²] (per cell).
        u: Wind speed [m/s] (per cell).
        tair_max: Maximum air temperature [K] (per cell).
        tair_min: Minimum air temperature [K] (per cell).
        qair: Relative humidity [0-1] (per cell).
        ch: Turbulent exchange coefficient [-] (per cell).
        alb: Surface albedo [-] (per cell).
        kaps: Soil thermal conductivity [W/m/K].
        nis: Soil thermal diffusivity [m²/s].
        tcost: Deep-soil constant temperature [K].
        pair: Air pressure [mb].
        ctim_s: Current-time seconds of day [s] (start of step).
        ftim_s: Future-time seconds of day [s] (end of step = ctim_s + dt).
        hrise_s: Sunrise time [s of day].
        hset_s: Sunset time [s of day].
        etrsuetp: Water-to-energy ET ratio [-]. ``1.0`` for the preliminary
            saturated pass; actual ``ETr/ETp`` for re-entry.
        dt: Model timestep [s].
        reentry: If True, skip the pre-sunrise night sub-period (physics is
            independent of ``etrsuetp``) and start day computation from
            ``td_rise`` rather than the current ``td``.

    Returns:
        ts_new (NDArray[np.float64]): Updated surface temperature [K].
        td_new (NDArray[np.float64]): Updated deep temperature [K].
        etp (NDArray[np.float64]): Potential evapotranspiration over the timestep [m].
        td_rise_new (NDArray[np.float64]): Updated ``td`` at sunrise (for later re-entry).
    """
    omega = 2.0 * np.pi / 86400.0

    # Decompose daily radiation into sinusoidal amplitude + constant.
    if dt > 86400.0 - 1.0:
        a_rad, c_rad = diurnal_radiation_cycle(rs, hrise_s, hset_s, "average")
    else:
        a_rad, c_rad = diurnal_radiation_cycle(rs, ctim_s, hset_s, "instant")

    # Air-temperature decomposition — constant over the day.
    a_tem = (tair_max - tair_min) / 2.0
    c_tem = (tair_max + tair_min) / 2.0  # tair_lin

    ts_out = ts.copy()
    td_out = td.copy()
    etp = np.zeros_like(ts)
    td_rise_out = td_rise.copy()

    def _night_call(
        td_in: NDArray[np.float64], t0: float, t1: float, kwargs: dict
    ) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
        return energy_balance_1l(
            omega,
            a_tem,
            0.0,
            _P_TEM_BASE + omega * t0,
            _P_RAD_BASE + omega * t0,
            c_tem,
            0.0,
            td_in,
            c_tem,
            0.0,
            pair,
            qair,
            t1 - t0,
            t1 - t0,
            ch,
            alb,
            kaps,
            nis,
            tcost,
            kwargs["etrsuetp"],
        )

    def _day_call(
        td_in: NDArray[np.float64], t0: float, t1: float, evap_step: float, kwargs: dict
    ) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
        return energy_balance_1l(
            omega,
            a_tem,
            a_rad,
            _P_TEM_BASE + omega * t0,
            _P_RAD_BASE + omega * t0,
            c_tem,
            c_rad,
            td_in,
            c_tem,
            u,
            pair,
            qair,
            t1 - t0,
            evap_step,
            ch,
            alb,
            kaps,
            nis,
            tcost,
            kwargs["etrsuetp"],
        )

    kw = {"etrsuetp": etrsuetp}

    if ctim_s < hrise_s:
        if ftim_s > hrise_s:
            # ctim < hrise < ftim  (step crosses sunrise)
            if not reentry:
                # Pre-sunrise night sub-period (independent of etrsuetp)
                ts_out, td_out, _ = _night_call(td_out, ctim_s, hrise_s, kw)
                td_rise_out = td_out.copy()
                day_start_td = td_out
            else:
                day_start_td = td_rise.copy()

            if ftim_s > hset_s:
                # ctim < hrise < hset < ftim  (full daylight inside step)
                ts_out, td_out, etp = _day_call(day_start_td, hrise_s, hset_s, _EVAP_SUBSTEP, kw)
                ts_out, td_out, _ = _night_call(td_out, hset_s, ftim_s, kw)
            else:
                # ctim < hrise < ftim < hset
                evap_step = min(ftim_s - hrise_s, _EVAP_SUBSTEP)
                ts_out, td_out, etp = _day_call(day_start_td, hrise_s, ftim_s, evap_step, kw)
        else:
            # ctim < ftim < hrise  (night only before sunrise)
            if reentry:
                # No radiation so etrsuetp has no impact — skip.
                return ts_out, td_out, etp, td_rise_out
            ts_out, td_out, _ = _night_call(td_out, ctim_s, ftim_s, kw)
    elif ctim_s < hset_s:
        if not reentry:
            td_rise_out = td_out.copy()
            day_start_td = td_out
        else:
            day_start_td = td_rise.copy()

        if ftim_s > hset_s:
            # hrise < ctim < hset < ftim  (step crosses sunset)
            evap_step = min(hset_s - ctim_s, _EVAP_SUBSTEP)
            ts_out, td_out, etp = _day_call(day_start_td, ctim_s, hset_s, evap_step, kw)
            ts_out, td_out, _ = _night_call(td_out, hset_s, ftim_s, kw)
        else:
            # hrise < ctim < ftim < hset  (daylight only)
            evap_step = min(ftim_s - ctim_s, _EVAP_SUBSTEP)
            ts_out, td_out, etp = _day_call(day_start_td, ctim_s, ftim_s, evap_step, kw)
    else:
        # hset < ctim < ftim  (night only after sunset)
        if reentry:
            return ts_out, td_out, etp, td_rise_out
        ts_out, td_out, _ = _night_call(td_out, ctim_s, ftim_s, kw)

    return ts_out, td_out, etp, td_rise_out

`mobidic.core.energy_balance.energy_balance_1l(ff, a_tem, a_rad, p_tem, p_rad, c_tem, c_rad, td_ini, tm, u, pair, hair, t_end, step, ch, alb, kaps, nis, tcost, etrsuetp)`

Analytical Fourier 1-layer surface energy balance over a sub-period.

Solves the constant + sinusoidal parts of the linearised surface energy balance analytically and integrates the evaporation flux via trapezoidal integration inside [0, t_end] with step step. Numba-accelerated.

Parameters:

Name	Type	Description	Default
`ff`	`float`	Diurnal angular frequency [rad/s] (typically `2*pi/86400`).	required
`a_tem`	`NDArray[float64]`	Air-temperature amplitude [K] (per cell).	required
`a_rad`	`NDArray[float64] \| float`	Radiation amplitude [W/m²] (per cell or scalar 0 for night).	required
`p_tem`	`float`	Air-temperature phase [rad] (scalar, already time-shifted).	required
`p_rad`	`float`	Radiation phase [rad] (scalar, already time-shifted).	required
`c_tem`	`NDArray[float64]`	Constant part of air temperature [K] (per cell).	required
`c_rad`	`NDArray[float64] \| float`	Constant part of radiation [W/m²] (per cell or 0 for night).	required
`td_ini`	`NDArray[float64]`	Initial deep-soil temperature [K] (per cell).	required
`tm`	`NDArray[float64]`	Mean air temperature used for linearisation [K] (per cell).	required
`u`	`NDArray[float64] \| float`	Wind speed [m/s] (per cell, or 0 for night).	required
`pair`	`float`	Air pressure [mb] (scalar).	required
`hair`	`NDArray[float64]`	Relative humidity [0-1] (per cell).	required
`t_end`	`float`	Integration length of this sub-period [s].	required
`step`	`float`	Integration step for evaporation [s].	required
`ch`	`NDArray[float64]`	Turbulent exchange coefficient for heat [-] (per cell).	required
`alb`	`NDArray[float64]`	Surface albedo [-] (per cell).	required
`kaps`	`float`	Soil thermal conductivity [W/m/K] (scalar).	required
`nis`	`float`	Soil thermal diffusivity [m²/s] (scalar).	required
`tcost`	`float`	Deep ground (constant) temperature [K] (scalar).	required
`etrsuetp`	`float \| NDArray[float64]`	Water-limited to energy-limited ET ratio [-] (scalar or per cell).	required

Returns:

Name	Type	Description
`ts`	`NDArray[float64]`	Surface temperature at the end of the sub-period [K].
`td`	`NDArray[float64]`	Deep-soil temperature at the end of the sub-period [K].
`evp`	`NDArray[float64]`	Evaporation over the sub-period [m].

Source code in mobidic/core/energy_balance.py

def energy_balance_1l(
    ff: float,
    a_tem: NDArray[np.float64],
    a_rad: NDArray[np.float64] | float,
    p_tem: float,
    p_rad: float,
    c_tem: NDArray[np.float64],
    c_rad: NDArray[np.float64] | float,
    td_ini: NDArray[np.float64],
    tm: NDArray[np.float64],
    u: NDArray[np.float64] | float,
    pair: float,
    hair: NDArray[np.float64],
    t_end: float,
    step: float,
    ch: NDArray[np.float64],
    alb: NDArray[np.float64],
    kaps: float,
    nis: float,
    tcost: float,
    etrsuetp: float | NDArray[np.float64],
) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
    """Analytical Fourier 1-layer surface energy balance over a sub-period.

    Solves the constant + sinusoidal parts of the linearised surface energy
    balance analytically and integrates the evaporation flux via trapezoidal
    integration inside ``[0, t_end]`` with step ``step``. Numba-accelerated.

    Args:
        ff: Diurnal angular frequency [rad/s] (typically ``2*pi/86400``).
        a_tem: Air-temperature amplitude [K] (per cell).
        a_rad: Radiation amplitude [W/m²] (per cell or scalar 0 for night).
        p_tem: Air-temperature phase [rad] (scalar, already time-shifted).
        p_rad: Radiation phase [rad] (scalar, already time-shifted).
        c_tem: Constant part of air temperature [K] (per cell).
        c_rad: Constant part of radiation [W/m²] (per cell or 0 for night).
        td_ini: Initial deep-soil temperature [K] (per cell).
        tm: Mean air temperature used for linearisation [K] (per cell).
        u: Wind speed [m/s] (per cell, or 0 for night).
        pair: Air pressure [mb] (scalar).
        hair: Relative humidity [0-1] (per cell).
        t_end: Integration length of this sub-period [s].
        step: Integration step for evaporation [s].
        ch: Turbulent exchange coefficient for heat [-] (per cell).
        alb: Surface albedo [-] (per cell).
        kaps: Soil thermal conductivity [W/m/K] (scalar).
        nis: Soil thermal diffusivity [m²/s] (scalar).
        tcost: Deep ground (constant) temperature [K] (scalar).
        etrsuetp: Water-limited to energy-limited ET ratio [-] (scalar or per cell).

    Returns:
        ts (NDArray[np.float64]): Surface temperature at the end of the sub-period [K].
        td (NDArray[np.float64]): Deep-soil temperature at the end of the sub-period [K].
        evp (NDArray[np.float64]): Evaporation over the sub-period [m].
    """
    tm_arr = np.ascontiguousarray(tm, dtype=np.float64)
    n = tm_arr.shape[0]

    a_tem_arr = _broadcast_to_n(a_tem, n)
    a_rad_arr = _broadcast_to_n(a_rad, n)
    c_tem_arr = _broadcast_to_n(c_tem, n)
    c_rad_arr = _broadcast_to_n(c_rad, n)
    td_ini_arr = _broadcast_to_n(td_ini, n)
    u_arr = _broadcast_to_n(u, n)
    hair_arr = _broadcast_to_n(hair, n)
    ch_arr = _broadcast_to_n(ch, n)
    alb_arr = _broadcast_to_n(alb, n)
    etr_arr = _broadcast_to_n(etrsuetp, n)

    if step <= 0:
        tt_values = np.zeros(1, dtype=np.float64)
        step_safe = 0.0
    else:
        n_pts = int(np.floor(t_end / step + 1e-9)) + 1
        tt_values = np.arange(n_pts, dtype=np.float64) * step
        step_safe = float(step)

    return _energy_balance_1l_kernel(
        float(ff),
        a_tem_arr,
        a_rad_arr,
        float(p_tem),
        float(p_rad),
        c_tem_arr,
        c_rad_arr,
        td_ini_arr,
        tm_arr,
        u_arr,
        float(pair),
        hair_arr,
        step_safe,
        ch_arr,
        alb_arr,
        float(kaps),
        float(nis),
        float(tcost),
        etr_arr,
        tt_values,
    )

`mobidic.core.energy_balance.diurnal_radiation_cycle(rs_avg, t_sunrise, t_sunset, mode='average')`

Decompose daily radiation into sinusoidal amplitude and constant.

The diurnal radiation cycle is represented as A * sin(omega * t + phi) + C where omega = 2*pi/86400 and the phase phi = -pi/2. This function returns (A, C).

Parameters:

Name	Type	Description	Default
`rs_avg`	`float \| NDArray[float64]`	Incoming shortwave radiation [W/m²]. Scalar or per-cell array.	required
`t_sunrise`	`float`	Sunrise time [s] (seconds of day).	required
`t_sunset`	`float`	Sunset time [s] (seconds of day).	required
`mode`	`str`	“average” if `rs_avg` is the day-average radiation, “instant” if it is an instantaneous value at `t_sunrise`.	`'average'`

Returns:

Type	Description
`tuple[NDArray[float64], NDArray[float64]]`	tuple[NDArray[np.float64], NDArray[np.float64]]: (amplitude, constant) arrays with the same shape as `rs_avg`.

Source code in mobidic/core/energy_balance.py

def diurnal_radiation_cycle(
    rs_avg: float | NDArray[np.float64],
    t_sunrise: float,
    t_sunset: float,
    mode: str = "average",
) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
    """Decompose daily radiation into sinusoidal amplitude and constant.

    The diurnal radiation cycle is represented as ``A * sin(omega * t + phi) + C``
    where ``omega = 2*pi/86400`` and the phase ``phi = -pi/2``. This function
    returns ``(A, C)``.

    Args:
        rs_avg: Incoming shortwave radiation [W/m²]. Scalar or per-cell array.
        t_sunrise: Sunrise time [s] (seconds of day).
        t_sunset: Sunset time [s] (seconds of day).
        mode: "average" if ``rs_avg`` is the day-average radiation,
            "instant" if it is an instantaneous value at ``t_sunrise``.

    Returns:
        tuple[NDArray[np.float64], NDArray[np.float64]]: (amplitude, constant) arrays with the same shape as ``rs_avg``.
    """
    w = 2.0 * np.pi / 86400.0
    p_r = -np.pi / 2.0

    if mode == "average":
        dd = (np.cos(w * t_sunrise + p_r) - np.cos(w * t_sunset + p_r)) / w - np.sin(w * t_sunset + p_r) * (
            t_sunset - t_sunrise
        )
        amplitude = rs_avg * (t_sunset - t_sunrise) / dd
    elif mode == "instant":
        dd = np.sin(w * t_sunrise + p_r) - np.sin(w * t_sunset + p_r)
        amplitude = rs_avg / dd
    else:
        raise ValueError(f"Invalid mode '{mode}', expected 'average' or 'instant'")

    constant = -amplitude * np.sin(w * t_sunset + p_r)
    return amplitude, constant

`mobidic.core.energy_balance.solar_hours(lat, lon, jday)`

Compute sunrise and sunset hours via bisection on solar elevation.

Parameters:

Name	Type	Description	Default
`lat`	`float`	Latitude [deg North].	required
`lon`	`float`	Longitude [deg East].	required
`jday`	`int`	Day of year (1-365).	required

Returns:

Type	Description
`tuple[float, float]`	tuple[float, float]: Sunrise (hrise) and sunset (hset) in decimal hours.

Source code in mobidic/core/energy_balance.py

def solar_hours(lat: float, lon: float, jday: int) -> tuple[float, float]:
    """Compute sunrise and sunset hours via bisection on solar elevation.

    Args:
        lat: Latitude [deg North].
        lon: Longitude [deg East].
        jday: Day of year (1-365).

    Returns:
        tuple[float, float]: Sunrise (hrise) and sunset (hset) in decimal hours.
    """
    # Sunrise: solar elevation crosses zero in the morning (0-12)
    hpre, hpos = 0.0, 12.0
    hrise = 6.0
    while (hpos - hpre) > (1.0 / 60.0):
        hrise = 0.5 * (hpos + hpre)
        _, el = solar_position(hrise, jday, lat, lon)
        if el > 0:
            hpos = hrise
        else:
            hpre = hrise

    # Sunset: solar elevation crosses zero in the afternoon (12-24)
    hpre, hpos = 12.0, 24.0
    hset = 18.0
    while (hpos - hpre) > (1.0 / 60.0):
        hset = 0.5 * (hpos + hpre)
        _, el = solar_position(hset, jday, lat, lon)
        if el > 0:
            hpre = hset
        else:
            hpos = hset

    return hrise, hset

`mobidic.core.energy_balance.solar_position(hour, day, lat, lon)`

Compute solar azimuth and elevation.

Parameters:

Name	Type	Description	Default
`hour`	`float`	Local time in decimal hours (0-24).	required
`day`	`int`	Day of year (Julian day, 1-365).	required
`lat`	`float`	Latitude [deg North].	required
`lon`	`float`	Longitude [deg East].	required

Returns:

Type	Description
`tuple[float, float]`	tuple[float, float]: Solar (azimuth, elevation) [deg].

Source code in mobidic/core/energy_balance.py

def solar_position(hour: float, day: int, lat: float, lon: float) -> tuple[float, float]:
    """Compute solar azimuth and elevation.

    Args:
        hour: Local time in decimal hours (0-24).
        day: Day of year (Julian day, 1-365).
        lat: Latitude [deg North].
        lon: Longitude [deg East].

    Returns:
        tuple[float, float]: Solar (azimuth, elevation) [deg].
    """
    delta = 23.45 * np.pi / 180.0 * np.cos(2.0 * np.pi / 365.0 * (172 - (day % 365)))
    dt1 = round(lon / 15.0) - lon / 15.0
    az = (hour + 12.0 - dt1) * 15.0
    az = az % 360.0
    sinalp = np.sin(delta) * np.sin(lat * np.pi / 180.0) + np.cos(delta) * np.cos(lat * np.pi / 180.0) * np.cos(
        az * np.pi / 180.0
    )
    el = 180.0 / np.pi * np.arcsin(sinalp)
    az = (180.0 + az) % 360.0
    return az, el

`mobidic.core.energy_balance.saturation_specific_humidity(T, P, dT=0.0)`

Saturation specific humidity (Magnus formula).

Parameters:

Name	Type	Description	Default
`T`	`NDArray[float64]`	Temperature [K].	required
`P`	`float`	Pressure [mb] (scalar).	required
`dT`	`float \| NDArray[float64]`	Dew-point depression [K] (default 0).	`0.0`

Returns:

Type	Description
`NDArray[float64]`	NDArray[np.float64]: Specific humidity [kg/kg].

Source code in mobidic/core/energy_balance.py

def saturation_specific_humidity(
    T: NDArray[np.float64],
    P: float,
    dT: float | NDArray[np.float64] = 0.0,
) -> NDArray[np.float64]:
    """Saturation specific humidity (Magnus formula).

    Args:
        T: Temperature [K].
        P: Pressure [mb] (scalar).
        dT: Dew-point depression [K] (default 0).

    Returns:
        NDArray[np.float64]: Specific humidity [kg/kg].
    """
    ep = 0.622
    Tc = T - 273.15 - dT
    es = 6.112 * np.exp(17.67 * Tc / (Tc + 243.5))
    # When dT is non-zero, apply the correction. dT=0 is the common case.
    if np.any(np.asarray(dT) != 0):
        es = es - (0.00066 * (1.0 + 0.00115 * Tc)) * P * dT
    return ep * es / (P - (1.0 - ep) * es)