Version: dev

SurfaceEnergyFlux

Net surface heat flux and its derivative with respect to surface temperature for the urban canopy thermal column. Migrated from UrbanCanopy.jl@9200edd4 src/roof_wall_road_snow_temperatures.jl:347 (`SurfaceEnergyFlux` @component). Implements Oleson et al. (2010, CLM Tech Note Ch. 4) Eqs. 4.26, 4.28, and 4.29: the surface heat flux h is the algebraic sum of absorbed solar radiation, net longwave radiation, sensible heat, latent heat, and waste/AC heat (Eq. 4.26); the longwave-derivative diagnostic dL_g_dT = 4*epsilon_g*sigma_SB*T_g_n^3 (Eq. 4.29) is the Stefan-Boltzmann temperature derivative of the upward longwave emission at the current surface temperature; the total flux derivative dh_dT = -dL_g_dT - dH_g_dT - dlambdaE_g_dT (Eq. 4.28) drives the linearized Crank-Nicholson scheme of Eq. 4.19. Parameterized as observed-only: S_g, L_g, H_g, lambdaE_g, H_wasteheat_g, H_aircond_g, epsilon_g, T_g_n, dH_g_dT, dlambdaE_g_dT are inputs (parameters); h, dL_g_dT, dh_dT are the algebraic state outputs.

.esm schema
1.0.0
Source
components/urban_canopy/temps/surface_energy_flux.esm

Description

Net surface heat flux and its derivative with respect to surface temperature for the urban canopy thermal column. Migrated from UrbanCanopy.jl@9200edd4 src/roof_wall_road_snow_temperatures.jl:347 (SurfaceEnergyFlux @component). Implements Oleson et al. (2010, CLM Tech Note Ch. 4) Eqs. 4.26, 4.28, and 4.29: the surface heat flux h is the algebraic sum of absorbed solar radiation, net longwave radiation, sensible heat, latent heat, and waste/AC heat (Eq. 4.26); the longwave-derivative diagnostic dL_g_dT = 4epsilon_gsigma_SB*T_g_n^3 (Eq. 4.29) is the Stefan-Boltzmann temperature derivative of the upward longwave emission at the current surface temperature; the total flux derivative dh_dT = -dL_g_dT - dH_g_dT - dlambdaE_g_dT (Eq. 4.28) drives the linearized Crank-Nicholson scheme of Eq. 4.19. Parameterized as observed-only: S_g, L_g, H_g, lambdaE_g, H_wasteheat_g, H_aircond_g, epsilon_g, T_g_n, dH_g_dT, dlambdaE_g_dT are inputs (parameters); h, dL_g_dT, dh_dT are the algebraic state outputs.

Reference

_source_ref: UrbanCanopy.jl@9200edd4f15609f20339d3a7d104de8064544310:src/roof_wall_road_snow_temperatures.jl:347 Upstream SurfaceEnergyFlux is an MTK @component with ten @parameters (S_g, L_g, H_g, λE_g, H_wasteheat_g, H_aircond_g, ε_g, T_g_n, dH_g_dT, dλE_g_dT), two @constants (σ_SB = 5.67e-8 W/(m^2*K^4) and one_K = 1.0 K), three time-indexed @variables (h, dh_dT, dL_g_dT), and three algebraic equations (no derivatives). Identical equation structure preserved here. Upstream canonical test (test/roof_wall_road_snow_temperatures_test.jl:282) drives the three equations through mtkcompile + ODEProblem + solve and asserts the three outputs at fixed inputs; that fixture is reproduced verbatim here as the upstream_equation_verification test.

The upstream component declares two @constants, both encoded here as parameters with the upstream defaults: sigma_SB = 5.67e-8 W/(m^2K^4) (Stefan-Boltzmann from Oleson Table 1.4) and one_K = 1.0 K (unit-balancer for the cube of T_g_n inside Eq. 4.29). Greek σ_SB, ε_g, λE_g, and dλE_g_dT are renamed sigma_SB, epsilon_g, lambdaE_g, and dlambdaE_g_dT for ASCII parameter names per the rig convention (matches the heating_cooling_flux.esm sibling’s encoding of upstream @constants α_CN → alpha_CN; snow_density.esm uses rho_).

Eq. 4.29 unit-balancing: the upstream form 4 * ε_g * σ_SB * (T_g_n / one_K)^3 * one_K^3 is preserved verbatim. The bare T_g_n^3 would not work in MTK’s symbolic algebra because the power operator demands a dimensionless exponent (here the literal 3) act on a dimensionless base; (T_g_n / one_K)^3 is dimensionless K/K cubed = dimensionless, and the trailing * one_K^3 restores the K^3 needed to balance the W/(m^2K^4) coefficient of 4*ε_g*σ_SB against the W/(m^2K) units of the LHS dL_g_dT. The runner observes the cube-of-cube structure as a chain of multiplies and one ^ op against the dimensionless quotient. The arithmetic value is identical to the unfused 4*ε_g*σ_SB*T_g_n^3 because one_K = 1.0 exactly, but the form is kept to mirror the upstream unit contract for any future ESS unit checker.

Eq. 4.26 sign convention: per the upstream description, h is the net heat flux INTO the surface (positive = energy gain). S_g is the absorbed solar radiation (positive = into the surface), L_g is the net longwave radiation in the convention where positive = net loss from the surface (consistent with the subtraction in the equation; the upstream description says downward positive but the equation form h = S - L - ... is loss-positive), H_g and lambdaE_g are upward sensible/latent fluxes (positive = loss from the surface), and the waste-heat and AC terms are additional energy inputs to the surface (positive = into the surface). The migrated ESM preserves the equation form exactly; the test negative_net_longwave_heat_gain pins the sign convention by exercising L_g < 0 (net longwave gain at the surface, h increases by |L_g|).

Observability classification: this component is functionally observed-only — the three output variables are each given by a single algebraic equation in terms of the parameters (the (t) time index has no role since there are no derivatives). Test scenarios use time_span: [0, 1] and sample at t=0 to mirror the convention used by sibling observed-only ESM files in this rig (heating_cooling_flux.esm, phase_change_energy.esm, total_phase_change_energy.esm).

Tolerance rationale (model-level rel=1e-9): observed-only algebraic system (no ODE integration); Eqs. 4.26 and 4.28 are additive at IEEE-754 ULP precision (~1e-16) while Eq. 4.29 involves one ^3 op with similar ULP-level roundoff. 1e-9 leaves comfortable headroom and matches sibling temps files (heating_cooling_flux.esm uses 1e-9, phase_change_energy.esm uses 1e-9). Tighter abs: 1e-12 is asserted on outputs that are exactly zero by construction (e.g. epsilon_g = 0 → dL_g_dT = 0; all-zero flux inputs → h = 0) so the assertion does not silently accept a stuck-at-1e-3 regression.

Upstream tolerance reproduction: upstream canonical test uses rtol = 1.0e-10 on h (purely additive) and rtol = 1.0e-6 on dL_g_dT and dh_dT (both involve the cube). Re-asserted at the same tolerances via per-assertion tolerance blocks on the upstream-canonical test, tighter than the model-level default for h._

Parameters

NameUnitsDefaultDescription
sigma_SBW/(m^2*K^4)5.67e-08Stefan-Boltzmann constant (Oleson Table 1.4) — upstream @constants σ_SB; the coefficient of the T^3 derivative of the Stefan-Boltzmann law for upward longwave emission used in Eq. 4.29
one_KK1.0Unit-temperature reference (1 K) — upstream @constants one_K; used in (T_g_n / one_K)^3 * one_K^3 to make the cube of T_g_n unit-consistent with MTK’s symbolic unit algebra (dimensionless base for the exponent, then restored K^3 for the outer dimension). Value is exactly 1.0 by construction; numerically a no-op.
S_gW/m^20.0Absorbed solar radiation at the surface (Oleson Eq. 4.26 first term) — net shortwave absorbed by the surface, positive into the surface
L_gW/m^20.0Net longwave radiation at the surface (Oleson Eq. 4.26 second term) — net upward longwave emission convention (positive = loss from the surface); subtracted in h
H_gW/m^20.0Sensible heat flux at the surface (Oleson Eq. 4.26 third term) — turbulent sensible heat loss from the surface to the atmosphere, positive = loss; subtracted in h
lambdaE_gW/m^20.0Latent heat flux at the surface (Oleson Eq. 4.26 fourth term) — upstream @parameters λE_g; turbulent latent heat loss from the surface to the atmosphere, positive = loss; subtracted in h
H_wasteheat_gW/m^20.0Waste heat flux applied to this surface (Oleson Eq. 4.26 fifth term) — anthropogenic waste-heat input from buildings, positive into the surface
H_aircond_gW/m^20.0Air-conditioning heat flux applied to this surface (Oleson Eq. 4.26 sixth term) — anthropogenic AC-driven heat rejection input, positive into the surface
epsilon_g10.95Surface emissivity (dimensionless) — upstream @parameters ε_g; coefficient on the Stefan-Boltzmann upward emission in Eq. 4.29. Typical urban surfaces 0.85-0.98; default 0.95 matches upstream test fixture
T_g_nK300.0Surface temperature at time n (Oleson Eq. 4.29) — input from the column-temperature solution at the start of the implicit-explicit step; cubed in the dL_g_dT derivative diagnostic
dH_g_dTW/(m^2*K)0.0Partial derivative of sensible heat flux with respect to T_g (Oleson Eq. 4.28) — supplied by an upstream turbulence/aerodynamic-resistance subsystem; enters dh_dT with a negative sign
dlambdaE_g_dTW/(m^2*K)0.0Partial derivative of latent heat flux with respect to T_g (Oleson Eq. 4.28) — upstream @parameters dλE_g_dT; supplied by the surface-flux subsystem (depends on saturation-vapor-pressure derivative); enters dh_dT with a negative sign

Observed

NameUnitsDefaultDescription
hW/m^2Net heat flux into the surface (Oleson Eq. 4.26): h = S_g - L_g - H_g - lambdaE_g + H_wasteheat_g + H_aircond_g; positive = net energy gain by the surface
dL_g_dTW/(m^2*K)Partial derivative of net longwave emission with respect to T_g (Oleson Eq. 4.29): 4epsilon_gsigma_SB*T_g_n^3; always non-negative because epsilon_g, sigma_SB, and T_g_n^3 are non-negative; characterizes the radiative-feedback strength of the surface temperature solver
dh_dTW/(m^2*K)Partial derivative of total surface heat flux with respect to T_g (Oleson Eq. 4.28): -dL_g_dT - dH_g_dT - dlambdaE_g_dT; drives the linearized Crank-Nicholson surface-temperature update (Eq. 4.19) — typically negative because warming the surface increases outgoing radiative and turbulent losses

Observed expressions

$$ h = S_{g} - L_{g} - H_{g} - lambdaE_{g} + H_{wasteheat\,g} + H_{aircond\,g} $$

Net heat flux into the surface (Oleson Eq. 4.26): h = S_g - L_g - H_g - lambdaE_g + H_wasteheat_g + H_aircond_g; positive = net energy gain by the surface

$$ dL_{g\,dT} = 4.0 \cdot epsilon_{g} \cdot sigma_{SB} \cdot \frac{T_{g\,n}}{one_{K}}^{3.0} \cdot one_{K}^{3.0} $$

Partial derivative of net longwave emission with respect to T_g (Oleson Eq. 4.29): 4epsilon_gsigma_SB*T_g_n^3; always non-negative because epsilon_g, sigma_SB, and T_g_n^3 are non-negative; characterizes the radiative-feedback strength of the surface temperature solver

$$ dh_{dT} = 0.0 - dL_{g\,dT} - dH_{g\,dT} - dlambdaE_{g\,dT} $$

Partial derivative of total surface heat flux with respect to T_g (Oleson Eq. 4.28): -dL_g_dT - dH_g_dT - dlambdaE_g_dT; drives the linearized Crank-Nicholson surface-temperature update (Eq. 4.19) — typically negative because warming the surface increases outgoing radiative and turbulent losses

Analyses

dL_g_dT_vs_T_g_n_cube_scaling

Sweep T_g_n across a realistic urban-surface temperature range (240-360 K) at fixed epsilon_g = 0.95 to surface the T^3 scaling of the longwave-derivative diagnostic (Eq. 4.29). dL_g_dT = 40.955.67e-8T_g_n^3 grows cubically with T_g_n, from 3.0 W/(m^2K) at 240 K to 10.0 W/(m^2K) at 360 K. Paired with the corresponding dh_dT curve (dh_dT = -dL_g_dT - constant turbulent terms): the constant offset of -7 W/(m^2K) (from dH_g_dT = 5 and dlambdaE_g_dT = 2) merely shifts the curve downward without changing its curvature.

Time span: 0.0 → 1.0

Swept over: T_g_n

dL_g_dT vs T_g_n: cubic curve through (0,0) with slope 4epsilon_gsigma_SB*T^2 — radiative feedback grows rapidly with surface temperature.

dh_dT vs T_g_n: -dL_g_dT - 7 — the same cube curve reflected across zero and shifted downward by the constant turbulent derivative offset.

dL_g_dT_heatmap_T_epsilon

Joint sweep over (T_g_n, epsilon_g) at fixed flux inputs, surfacing the bilinear-cubic structure of Eq. 4.29: dL_g_dT scales linearly with epsilon_g (vertical axis) and cubically with T_g_n (horizontal axis). Iso-dL_g_dT contours follow curves of constant epsilon_g * T_g_n^3 — at low T_g_n a wide range of emissivities give similar small feedback; at high T_g_n the contour spacing tightens, showing the strong sensitivity to T in summer urban heat-island regimes. The vertical strip at epsilon_g = 0 collapses to dL_g_dT = 0 regardless of T_g_n.

Time span: 0.0 → 1.0

Swept over: T_g_n, epsilon_g

dL_g_dT in (T_g_n, epsilon_g) space: bilinear-cubic surface with iso-contours along constant epsilon_g * T_g_n^3.

Raw .esm

Raw .esm JSON (this component)
{
  "reference": {
    "notes": "source_ref: UrbanCanopy.jl@9200edd4f15609f20339d3a7d104de8064544310:src/roof_wall_road_snow_temperatures.jl:347\nUpstream `SurfaceEnergyFlux` is an MTK `@component` with ten `@parameters` (S_g, L_g, H_g, λE_g, H_wasteheat_g, H_aircond_g, ε_g, T_g_n, dH_g_dT, dλE_g_dT), two `@constants` (σ_SB = 5.67e-8 W/(m^2*K^4) and one_K = 1.0 K), three time-indexed `@variables` (h, dh_dT, dL_g_dT), and three algebraic equations (no derivatives). Identical equation structure preserved here. Upstream canonical test (test/roof_wall_road_snow_temperatures_test.jl:282) drives the three equations through `mtkcompile` + `ODEProblem` + `solve` and asserts the three outputs at fixed inputs; that fixture is reproduced verbatim here as the `upstream_equation_verification` test.\n\nThe upstream component declares two `@constants`, both encoded here as parameters with the upstream defaults: sigma_SB = 5.67e-8 W/(m^2*K^4) (Stefan-Boltzmann from Oleson Table 1.4) and one_K = 1.0 K (unit-balancer for the cube of T_g_n inside Eq. 4.29). Greek `σ_SB`, `ε_g`, `λE_g`, and `dλE_g_dT` are renamed `sigma_SB`, `epsilon_g`, `lambdaE_g`, and `dlambdaE_g_dT` for ASCII parameter names per the rig convention (matches the heating_cooling_flux.esm sibling's encoding of upstream `@constants` α_CN → alpha_CN; snow_density.esm uses rho_*).\n\nEq. 4.29 unit-balancing: the upstream form `4 * ε_g * σ_SB * (T_g_n / one_K)^3 * one_K^3` is preserved verbatim. The bare `T_g_n^3` would not work in MTK's symbolic algebra because the power operator demands a dimensionless exponent (here the literal `3`) act on a dimensionless base; `(T_g_n / one_K)^3` is dimensionless K/K cubed = dimensionless, and the trailing `* one_K^3` restores the K^3 needed to balance the W/(m^2*K^4) coefficient of `4*ε_g*σ_SB` against the W/(m^2*K) units of the LHS `dL_g_dT`. The runner observes the cube-of-cube structure as a chain of multiplies and one `^` op against the dimensionless quotient. The arithmetic value is identical to the unfused `4*ε_g*σ_SB*T_g_n^3` because one_K = 1.0 exactly, but the form is kept to mirror the upstream unit contract for any future ESS unit checker.\n\nEq. 4.26 sign convention: per the upstream description, h is the net heat flux INTO the surface (positive = energy gain). S_g is the absorbed solar radiation (positive = into the surface), L_g is the net longwave radiation in the convention where positive = net loss from the surface (consistent with the subtraction in the equation; the upstream description says `downward positive` but the equation form `h = S - L - ...` is loss-positive), H_g and lambdaE_g are upward sensible/latent fluxes (positive = loss from the surface), and the waste-heat and AC terms are additional energy inputs to the surface (positive = into the surface). The migrated ESM preserves the equation form exactly; the test `negative_net_longwave_heat_gain` pins the sign convention by exercising L_g < 0 (net longwave gain at the surface, h increases by |L_g|).\n\nObservability classification: this component is functionally observed-only — the three output variables are each given by a single algebraic equation in terms of the parameters (the `(t)` time index has no role since there are no derivatives). Test scenarios use `time_span: [0, 1]` and sample at t=0 to mirror the convention used by sibling observed-only ESM files in this rig (heating_cooling_flux.esm, phase_change_energy.esm, total_phase_change_energy.esm).\n\nTolerance rationale (model-level rel=1e-9): observed-only algebraic system (no ODE integration); Eqs. 4.26 and 4.28 are additive at IEEE-754 ULP precision (~1e-16) while Eq. 4.29 involves one `^3` op with similar ULP-level roundoff. 1e-9 leaves comfortable headroom and matches sibling temps files (heating_cooling_flux.esm uses 1e-9, phase_change_energy.esm uses 1e-9). Tighter `abs: 1e-12` is asserted on outputs that are exactly zero by construction (e.g. `epsilon_g = 0` → dL_g_dT = 0; all-zero flux inputs → h = 0) so the assertion does not silently accept a stuck-at-1e-3 regression.\n\nUpstream tolerance reproduction: upstream canonical test uses `rtol = 1.0e-10` on h (purely additive) and `rtol = 1.0e-6` on dL_g_dT and dh_dT (both involve the cube). Re-asserted at the same tolerances via per-assertion `tolerance` blocks on the upstream-canonical test, tighter than the model-level default for h."
  },
  "tolerance": {
    "rel": 1e-09
  },
  "variables": {
    "sigma_SB": {
      "type": "parameter",
      "units": "W/(m^2*K^4)",
      "default": 5.67e-08,
      "description": "Stefan-Boltzmann constant (Oleson Table 1.4) — upstream `@constants` σ_SB; the coefficient of the T^3 derivative of the Stefan-Boltzmann law for upward longwave emission used in Eq. 4.29"
    },
    "one_K": {
      "type": "parameter",
      "units": "K",
      "default": 1.0,
      "description": "Unit-temperature reference (1 K) — upstream `@constants` one_K; used in `(T_g_n / one_K)^3 * one_K^3` to make the cube of T_g_n unit-consistent with MTK's symbolic unit algebra (dimensionless base for the exponent, then restored K^3 for the outer dimension). Value is exactly 1.0 by construction; numerically a no-op."
    },
    "S_g": {
      "type": "parameter",
      "units": "W/m^2",
      "default": 0.0,
      "description": "Absorbed solar radiation at the surface (Oleson Eq. 4.26 first term) — net shortwave absorbed by the surface, positive into the surface"
    },
    "L_g": {
      "type": "parameter",
      "units": "W/m^2",
      "default": 0.0,
      "description": "Net longwave radiation at the surface (Oleson Eq. 4.26 second term) — net upward longwave emission convention (positive = loss from the surface); subtracted in h"
    },
    "H_g": {
      "type": "parameter",
      "units": "W/m^2",
      "default": 0.0,
      "description": "Sensible heat flux at the surface (Oleson Eq. 4.26 third term) — turbulent sensible heat loss from the surface to the atmosphere, positive = loss; subtracted in h"
    },
    "lambdaE_g": {
      "type": "parameter",
      "units": "W/m^2",
      "default": 0.0,
      "description": "Latent heat flux at the surface (Oleson Eq. 4.26 fourth term) — upstream `@parameters` λE_g; turbulent latent heat loss from the surface to the atmosphere, positive = loss; subtracted in h"
    },
    "H_wasteheat_g": {
      "type": "parameter",
      "units": "W/m^2",
      "default": 0.0,
      "description": "Waste heat flux applied to this surface (Oleson Eq. 4.26 fifth term) — anthropogenic waste-heat input from buildings, positive into the surface"
    },
    "H_aircond_g": {
      "type": "parameter",
      "units": "W/m^2",
      "default": 0.0,
      "description": "Air-conditioning heat flux applied to this surface (Oleson Eq. 4.26 sixth term) — anthropogenic AC-driven heat rejection input, positive into the surface"
    },
    "epsilon_g": {
      "type": "parameter",
      "units": "1",
      "default": 0.95,
      "description": "Surface emissivity (dimensionless) — upstream `@parameters` ε_g; coefficient on the Stefan-Boltzmann upward emission in Eq. 4.29. Typical urban surfaces 0.85-0.98; default 0.95 matches upstream test fixture"
    },
    "T_g_n": {
      "type": "parameter",
      "units": "K",
      "default": 300.0,
      "description": "Surface temperature at time n (Oleson Eq. 4.29) — input from the column-temperature solution at the start of the implicit-explicit step; cubed in the dL_g_dT derivative diagnostic"
    },
    "dH_g_dT": {
      "type": "parameter",
      "units": "W/(m^2*K)",
      "default": 0.0,
      "description": "Partial derivative of sensible heat flux with respect to T_g (Oleson Eq. 4.28) — supplied by an upstream turbulence/aerodynamic-resistance subsystem; enters dh_dT with a negative sign"
    },
    "dlambdaE_g_dT": {
      "type": "parameter",
      "units": "W/(m^2*K)",
      "default": 0.0,
      "description": "Partial derivative of latent heat flux with respect to T_g (Oleson Eq. 4.28) — upstream `@parameters` dλE_g_dT; supplied by the surface-flux subsystem (depends on saturation-vapor-pressure derivative); enters dh_dT with a negative sign"
    },
    "h": {
      "type": "unknown",
      "units": "W/m^2",
      "description": "Net heat flux into the surface (Oleson Eq. 4.26): h = S_g - L_g - H_g - lambdaE_g + H_wasteheat_g + H_aircond_g; positive = net energy gain by the surface"
    },
    "dL_g_dT": {
      "type": "unknown",
      "units": "W/(m^2*K)",
      "description": "Partial derivative of net longwave emission with respect to T_g (Oleson Eq. 4.29): 4*epsilon_g*sigma_SB*T_g_n^3; always non-negative because epsilon_g, sigma_SB, and T_g_n^3 are non-negative; characterizes the radiative-feedback strength of the surface temperature solver"
    },
    "dh_dT": {
      "type": "unknown",
      "units": "W/(m^2*K)",
      "description": "Partial derivative of total surface heat flux with respect to T_g (Oleson Eq. 4.28): -dL_g_dT - dH_g_dT - dlambdaE_g_dT; drives the linearized Crank-Nicholson surface-temperature update (Eq. 4.19) — typically negative because warming the surface increases outgoing radiative and turbulent losses"
    }
  },
  "equations": [
    {
      "lhs": "h",
      "rhs": {
        "op": "+",
        "args": [
          {
            "op": "+",
            "args": [
              {
                "op": "-",
                "args": [
                  {
                    "op": "-",
                    "args": [
                      {
                        "op": "-",
                        "args": [
                          "S_g",
                          "L_g"
                        ]
                      },
                      "H_g"
                    ]
                  },
                  "lambdaE_g"
                ]
              },
              "H_wasteheat_g"
            ]
          },
          "H_aircond_g"
        ]
      }
    },
    {
      "lhs": "dL_g_dT",
      "rhs": {
        "op": "*",
        "args": [
          {
            "op": "*",
            "args": [
              {
                "op": "*",
                "args": [
                  {
                    "op": "*",
                    "args": [
                      4.0,
                      "epsilon_g"
                    ]
                  },
                  "sigma_SB"
                ]
              },
              {
                "op": "^",
                "args": [
                  {
                    "op": "/",
                    "args": [
                      "T_g_n",
                      "one_K"
                    ]
                  },
                  3.0
                ]
              }
            ]
          },
          {
            "op": "^",
            "args": [
              "one_K",
              3.0
            ]
          }
        ]
      }
    },
    {
      "lhs": "dh_dT",
      "rhs": {
        "op": "-",
        "args": [
          {
            "op": "-",
            "args": [
              {
                "op": "-",
                "args": [
                  0.0,
                  "dL_g_dT"
                ]
              },
              "dH_g_dT"
            ]
          },
          "dlambdaE_g_dT"
        ]
      }
    }
  ],
  "tests": [
    {
      "id": "upstream_equation_verification",
      "description": "Reproduces upstream `SurfaceEnergyFlux - Equation Verification` test (UrbanCanopy.jl test/roof_wall_road_snow_temperatures_test.jl:282). Canonical inputs: S_g = 200, L_g = 50, H_g = 30, lambdaE_g = 20, H_wasteheat_g = 10, H_aircond_g = 5, epsilon_g = 0.95, T_g_n = 300, dH_g_dT = 5, dlambdaE_g_dT = 2. Expected: h = 200 - 50 - 30 - 20 + 10 + 5 = 115 W/m^2 (Eq. 4.26); dL_g_dT = 4*0.95*5.67e-8*300^3 = 5.81742 W/(m^2*K) (Eq. 4.29); dh_dT = -5.81742 - 5 - 2 = -12.81742 W/(m^2*K) (Eq. 4.28). Tolerances mirror upstream: rtol = 1.0e-10 on h (purely additive) and rtol = 1.0e-6 on the two cube-based outputs.",
      "parameter_overrides": {
        "S_g": 200.0,
        "L_g": 50.0,
        "H_g": 30.0,
        "lambdaE_g": 20.0,
        "H_wasteheat_g": 10.0,
        "H_aircond_g": 5.0,
        "epsilon_g": 0.95,
        "T_g_n": 300.0,
        "dH_g_dT": 5.0,
        "dlambdaE_g_dT": 2.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "h",
          "time": 0.0,
          "expected": 115.0,
          "tolerance": {
            "rel": 1e-10
          }
        },
        {
          "variable": "dL_g_dT",
          "time": 0.0,
          "expected": 5.81742,
          "tolerance": {
            "rel": 1e-06
          }
        },
        {
          "variable": "dh_dT",
          "time": 0.0,
          "expected": -12.81742,
          "tolerance": {
            "rel": 1e-06
          }
        }
      ]
    },
    {
      "id": "all_zero_inputs",
      "description": "Degenerate boundary: every flux parameter and derivative is zero, T_g_n = 0 K, epsilon_g = 0. h = 0 (additive identity), dL_g_dT = 4 * 0 * sigma_SB * 0^3 = 0 exactly, dh_dT = -0 - 0 - 0 = 0. Anchors the all-zero result of the merged equations and exercises the runner against epsilon_g = 0 and T_g_n = 0 simultaneously (each independently zeros dL_g_dT — verifies neither raises an underflow / divide-by-zero error).",
      "parameter_overrides": {
        "S_g": 0.0,
        "L_g": 0.0,
        "H_g": 0.0,
        "lambdaE_g": 0.0,
        "H_wasteheat_g": 0.0,
        "H_aircond_g": 0.0,
        "epsilon_g": 0.0,
        "T_g_n": 0.0,
        "dH_g_dT": 0.0,
        "dlambdaE_g_dT": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "h",
          "time": 0.0,
          "expected": 0.0,
          "tolerance": {
            "abs": 1e-12
          }
        },
        {
          "variable": "dL_g_dT",
          "time": 0.0,
          "expected": 0.0,
          "tolerance": {
            "abs": 1e-12
          }
        },
        {
          "variable": "dh_dT",
          "time": 0.0,
          "expected": 0.0,
          "tolerance": {
            "abs": 1e-12
          }
        }
      ]
    },
    {
      "id": "zero_emissivity_no_radiative_feedback",
      "description": "Boundary case: epsilon_g = 0 (perfectly reflective, hypothetical limit). The longwave-derivative term vanishes: dL_g_dT = 4 * 0 * sigma_SB * T_g_n^3 * 1 = 0 exactly regardless of T_g_n. The radiative feedback in dh_dT therefore collapses to dh_dT = -dH_g_dT - dlambdaE_g_dT = -5 - 2 = -7 W/(m^2*K). h is unaffected by epsilon_g (it does not appear in Eq. 4.26). Isolates the epsilon_g pre-factor of Eq. 4.29.",
      "parameter_overrides": {
        "S_g": 200.0,
        "L_g": 50.0,
        "H_g": 30.0,
        "lambdaE_g": 20.0,
        "H_wasteheat_g": 10.0,
        "H_aircond_g": 5.0,
        "epsilon_g": 0.0,
        "T_g_n": 300.0,
        "dH_g_dT": 5.0,
        "dlambdaE_g_dT": 2.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "h",
          "time": 0.0,
          "expected": 115.0
        },
        {
          "variable": "dL_g_dT",
          "time": 0.0,
          "expected": 0.0,
          "tolerance": {
            "abs": 1e-12
          }
        },
        {
          "variable": "dh_dT",
          "time": 0.0,
          "expected": -7.0
        }
      ]
    },
    {
      "id": "blackbody_epsilon_one",
      "description": "Boundary case: epsilon_g = 1 (blackbody upper limit). dL_g_dT = 4 * 1 * 5.67e-8 * 300^3 = 6.1236 W/(m^2*K) — the largest possible radiative feedback at T_g_n = 300 K. dh_dT = -6.1236 - 5 - 2 = -13.1236 W/(m^2*K). Anchors the upper bound of the radiative-feedback magnitude and verifies that the equation is linear in epsilon_g (the ratio dL_g_dT(epsilon=1) / dL_g_dT(epsilon=0.95) = 1/0.95 ≈ 1.0526; both are reported, so a reviewer can cross-check).",
      "parameter_overrides": {
        "S_g": 200.0,
        "L_g": 50.0,
        "H_g": 30.0,
        "lambdaE_g": 20.0,
        "H_wasteheat_g": 10.0,
        "H_aircond_g": 5.0,
        "epsilon_g": 1.0,
        "T_g_n": 300.0,
        "dH_g_dT": 5.0,
        "dlambdaE_g_dT": 2.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "h",
          "time": 0.0,
          "expected": 115.0
        },
        {
          "variable": "dL_g_dT",
          "time": 0.0,
          "expected": 6.1236,
          "tolerance": {
            "rel": 1e-06
          }
        },
        {
          "variable": "dh_dT",
          "time": 0.0,
          "expected": -13.1236,
          "tolerance": {
            "rel": 1e-06
          }
        }
      ]
    },
    {
      "id": "zero_surface_temperature_no_radiation",
      "description": "Boundary case: T_g_n = 0 K (unphysical but mathematically valid in the equation). dL_g_dT = 4 * 0.95 * 5.67e-8 * 0^3 * 1 = 0 exactly — the cube zeroes out the entire radiative-derivative term regardless of epsilon_g. dh_dT = -0 - dH_g_dT - dlambdaE_g_dT = -5 - 2 = -7 W/(m^2*K). Exercises the cube power against zero and confirms IEEE-754 0^3 = 0 (no NaN/Inf branch).",
      "parameter_overrides": {
        "S_g": 200.0,
        "L_g": 50.0,
        "H_g": 30.0,
        "lambdaE_g": 20.0,
        "H_wasteheat_g": 10.0,
        "H_aircond_g": 5.0,
        "epsilon_g": 0.95,
        "T_g_n": 0.0,
        "dH_g_dT": 5.0,
        "dlambdaE_g_dT": 2.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "h",
          "time": 0.0,
          "expected": 115.0
        },
        {
          "variable": "dL_g_dT",
          "time": 0.0,
          "expected": 0.0,
          "tolerance": {
            "abs": 1e-12
          }
        },
        {
          "variable": "dh_dT",
          "time": 0.0,
          "expected": -7.0
        }
      ]
    },
    {
      "id": "cold_surface_lower_radiative_feedback",
      "description": "Realistic cold-surface case: T_g_n = 250 K (e.g. nighttime urban canyon in winter). dL_g_dT = 4 * 0.95 * 5.67e-8 * 250^3 = 3.3665625 W/(m^2*K) — substantially smaller than at 300 K because of the cube scaling (T^3 ratio = (250/300)^3 ≈ 0.579). dh_dT = -3.3665625 - 5 - 2 = -10.3665625 W/(m^2*K). Surfaces the cube-of-temperature scaling that dominates the seasonal variation of the radiative feedback term.",
      "parameter_overrides": {
        "S_g": 200.0,
        "L_g": 50.0,
        "H_g": 30.0,
        "lambdaE_g": 20.0,
        "H_wasteheat_g": 10.0,
        "H_aircond_g": 5.0,
        "epsilon_g": 0.95,
        "T_g_n": 250.0,
        "dH_g_dT": 5.0,
        "dlambdaE_g_dT": 2.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "h",
          "time": 0.0,
          "expected": 115.0
        },
        {
          "variable": "dL_g_dT",
          "time": 0.0,
          "expected": 3.3665625,
          "tolerance": {
            "rel": 1e-06
          }
        },
        {
          "variable": "dh_dT",
          "time": 0.0,
          "expected": -10.3665625,
          "tolerance": {
            "rel": 1e-06
          }
        }
      ]
    },
    {
      "id": "hot_surface_higher_radiative_feedback",
      "description": "Realistic hot-surface case: T_g_n = 350 K (e.g. midsummer sunlit asphalt). dL_g_dT = 4 * 0.95 * 5.67e-8 * 350^3 = 9.2378475 W/(m^2*K) — substantially larger than at 300 K. dh_dT = -9.2378475 - 5 - 2 = -16.2378475 W/(m^2*K). Pairs with `cold_surface_lower_radiative_feedback` to surface the monotone increase of |dh_dT| with T_g_n (warmer surfaces have stronger negative feedback, stabilizing the implicit Crank-Nicholson update).",
      "parameter_overrides": {
        "S_g": 200.0,
        "L_g": 50.0,
        "H_g": 30.0,
        "lambdaE_g": 20.0,
        "H_wasteheat_g": 10.0,
        "H_aircond_g": 5.0,
        "epsilon_g": 0.95,
        "T_g_n": 350.0,
        "dH_g_dT": 5.0,
        "dlambdaE_g_dT": 2.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "h",
          "time": 0.0,
          "expected": 115.0
        },
        {
          "variable": "dL_g_dT",
          "time": 0.0,
          "expected": 9.2378475,
          "tolerance": {
            "rel": 1e-06
          }
        },
        {
          "variable": "dh_dT",
          "time": 0.0,
          "expected": -16.2378475,
          "tolerance": {
            "rel": 1e-06
          }
        }
      ]
    },
    {
      "id": "wasteheat_and_aircond_isolated",
      "description": "Isolates the anthropogenic-heat terms: all natural fluxes zero (S_g = L_g = H_g = lambdaE_g = 0) with H_wasteheat_g = 50 W/m^2 and H_aircond_g = 25 W/m^2. h = 0 - 0 - 0 - 0 + 50 + 25 = 75 W/m^2. epsilon_g and T_g_n still drive dL_g_dT and dh_dT as before. Anchors the additive contribution of the urban anthropogenic terms in Eq. 4.26 — these distinguish the CLMU urban-canopy energy balance from the underlying CLM natural-surface energy balance.",
      "parameter_overrides": {
        "S_g": 0.0,
        "L_g": 0.0,
        "H_g": 0.0,
        "lambdaE_g": 0.0,
        "H_wasteheat_g": 50.0,
        "H_aircond_g": 25.0,
        "epsilon_g": 0.95,
        "T_g_n": 300.0,
        "dH_g_dT": 0.0,
        "dlambdaE_g_dT": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "h",
          "time": 0.0,
          "expected": 75.0
        },
        {
          "variable": "dL_g_dT",
          "time": 0.0,
          "expected": 5.81742,
          "tolerance": {
            "rel": 1e-06
          }
        },
        {
          "variable": "dh_dT",
          "time": 0.0,
          "expected": -5.81742,
          "tolerance": {
            "rel": 1e-06
          }
        }
      ]
    },
    {
      "id": "negative_net_longwave_heat_gain",
      "description": "Regime check: L_g = -30 W/m^2 (more incoming longwave than outgoing — e.g. cloudy night surface emitting less than the atmosphere emits downward). The subtraction `-L_g = -(-30) = +30` adds to h, giving h = 200 - (-30) - 30 - 20 + 10 + 5 = 195 W/m^2. Confirms the sign convention of Eq. 4.26: negative L_g represents net longwave gain at the surface and is correctly handled by the signed subtraction.",
      "parameter_overrides": {
        "S_g": 200.0,
        "L_g": -30.0,
        "H_g": 30.0,
        "lambdaE_g": 20.0,
        "H_wasteheat_g": 10.0,
        "H_aircond_g": 5.0,
        "epsilon_g": 0.95,
        "T_g_n": 300.0,
        "dH_g_dT": 5.0,
        "dlambdaE_g_dT": 2.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "h",
          "time": 0.0,
          "expected": 195.0
        },
        {
          "variable": "dL_g_dT",
          "time": 0.0,
          "expected": 5.81742,
          "tolerance": {
            "rel": 1e-06
          }
        }
      ]
    },
    {
      "id": "zero_turbulent_derivatives_only_radiative_feedback",
      "description": "dH_g_dT = dlambdaE_g_dT = 0 isolates the radiative term in Eq. 4.28: dh_dT = -dL_g_dT - 0 - 0 = -5.81742 W/(m^2*K) at canonical T_g_n = 300, epsilon_g = 0.95. Exercises the no-turbulence limit of the surface-temperature solver feedback and confirms additive structure of Eq. 4.28.",
      "parameter_overrides": {
        "S_g": 200.0,
        "L_g": 50.0,
        "H_g": 30.0,
        "lambdaE_g": 20.0,
        "H_wasteheat_g": 10.0,
        "H_aircond_g": 5.0,
        "epsilon_g": 0.95,
        "T_g_n": 300.0,
        "dH_g_dT": 0.0,
        "dlambdaE_g_dT": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "h",
          "time": 0.0,
          "expected": 115.0
        },
        {
          "variable": "dL_g_dT",
          "time": 0.0,
          "expected": 5.81742,
          "tolerance": {
            "rel": 1e-06
          }
        },
        {
          "variable": "dh_dT",
          "time": 0.0,
          "expected": -5.81742,
          "tolerance": {
            "rel": 1e-06
          }
        }
      ]
    }
  ],
  "analyses": [
    {
      "id": "dL_g_dT_vs_T_g_n_cube_scaling",
      "description": "Sweep T_g_n across a realistic urban-surface temperature range (240-360 K) at fixed epsilon_g = 0.95 to surface the T^3 scaling of the longwave-derivative diagnostic (Eq. 4.29). dL_g_dT = 4*0.95*5.67e-8*T_g_n^3 grows cubically with T_g_n, from 3.0 W/(m^2*K) at 240 K to 10.0 W/(m^2*K) at 360 K. Paired with the corresponding dh_dT curve (dh_dT = -dL_g_dT - constant turbulent terms): the constant offset of -7 W/(m^2*K) (from dH_g_dT = 5 and dlambdaE_g_dT = 2) merely shifts the curve downward without changing its curvature.",
      "initial_state": {},
      "parameters": {
        "S_g": 200.0,
        "L_g": 50.0,
        "H_g": 30.0,
        "lambdaE_g": 20.0,
        "H_wasteheat_g": 10.0,
        "H_aircond_g": 5.0,
        "epsilon_g": 0.95,
        "dH_g_dT": 5.0,
        "dlambdaE_g_dT": 2.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "parameter_sweep": {
        "type": "cartesian",
        "dimensions": [
          {
            "parameter": "T_g_n",
            "range": {
              "start": 240.0,
              "stop": 360.0,
              "count": 31,
              "scale": "linear"
            }
          }
        ]
      },
      "plots": [
        {
          "id": "dL_g_dT_curve",
          "type": "line",
          "description": "dL_g_dT vs T_g_n: cubic curve through (0,0) with slope 4*epsilon_g*sigma_SB*T^2 — radiative feedback grows rapidly with surface temperature.",
          "x": {
            "variable": "T_g_n",
            "label": "T_g_n (K)"
          },
          "y": {
            "variable": "dL_g_dT",
            "label": "dL_g_dT (W/(m^2*K))"
          }
        },
        {
          "id": "dh_dT_curve",
          "type": "line",
          "description": "dh_dT vs T_g_n: -dL_g_dT - 7 — the same cube curve reflected across zero and shifted downward by the constant turbulent derivative offset.",
          "x": {
            "variable": "T_g_n",
            "label": "T_g_n (K)"
          },
          "y": {
            "variable": "dh_dT",
            "label": "dh_dT (W/(m^2*K))"
          }
        }
      ]
    },
    {
      "id": "dL_g_dT_heatmap_T_epsilon",
      "description": "Joint sweep over (T_g_n, epsilon_g) at fixed flux inputs, surfacing the bilinear-cubic structure of Eq. 4.29: dL_g_dT scales linearly with epsilon_g (vertical axis) and cubically with T_g_n (horizontal axis). Iso-dL_g_dT contours follow curves of constant epsilon_g * T_g_n^3 — at low T_g_n a wide range of emissivities give similar small feedback; at high T_g_n the contour spacing tightens, showing the strong sensitivity to T in summer urban heat-island regimes. The vertical strip at epsilon_g = 0 collapses to dL_g_dT = 0 regardless of T_g_n.",
      "initial_state": {},
      "parameters": {
        "S_g": 200.0,
        "L_g": 50.0,
        "H_g": 30.0,
        "lambdaE_g": 20.0,
        "H_wasteheat_g": 10.0,
        "H_aircond_g": 5.0,
        "dH_g_dT": 5.0,
        "dlambdaE_g_dT": 2.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "parameter_sweep": {
        "type": "cartesian",
        "dimensions": [
          {
            "parameter": "T_g_n",
            "range": {
              "start": 240.0,
              "stop": 360.0,
              "count": 13,
              "scale": "linear"
            }
          },
          {
            "parameter": "epsilon_g",
            "range": {
              "start": 0.0,
              "stop": 1.0,
              "count": 11,
              "scale": "linear"
            }
          }
        ]
      },
      "plots": [
        {
          "id": "dL_g_dT_T_epsilon_heatmap",
          "type": "heatmap",
          "description": "dL_g_dT in (T_g_n, epsilon_g) space: bilinear-cubic surface with iso-contours along constant epsilon_g * T_g_n^3.",
          "x": {
            "variable": "T_g_n",
            "label": "T_g_n (K)"
          },
          "y": {
            "variable": "epsilon_g",
            "label": "epsilon_g (dimensionless)"
          },
          "value": {
            "variable": "dL_g_dT",
            "reduce": "final"
          }
        }
      ]
    }
  ]
}