PhaseChangeEnergy

{
  "reference": {
    "notes": "source_ref: UrbanCanopy.jl@9200edd4f15609f20339d3a7d104de8064544310:src/roof_wall_road_snow_temperatures.jl:496\nUpstream `PhaseChangeEnergy` is an MTK `@component` keyword-parameterized by `layer_type ∈ {:interior, :top}` (default `:interior`). It declares three `@constants` (T_f = 273.15 K, L_f = 3.337e5 J/kg, α_CN = 0.5), four shared `@parameters` (c_i, Δz_i, Δt, T_i_n), two shared bottom-flux `@parameters` (F_i_n, F_i_np1), one time-indexed `@variable` (H_excess), and one algebraic equation. The `:top` branch adds two `@parameters` (h_n, dh_dT) and uses Eq. 4.59 for i = snl+1; the `:interior` branch adds two `@parameters` (F_im1_n, F_im1_np1) and uses Eq. 4.59 for i > snl+1.\n\nVariant unification: rather than emitting two separate ESM files for the two upstream variants (no convention in this rig for a single tracker row mapped to multiple files), the two equation variants are merged into a single algebraic equation gated by a dimensionless switch parameter `is_top` (default 0.0 = interior layer). The combined form is\n\n    H_excess = α_CN*F_i_n + (1-α_CN)*F_i_np1\n               + is_top * (h_n + dh_dT*(T_f - T_i_n))\n               - (1 - is_top) * (α_CN*F_im1_n + (1-α_CN)*F_im1_np1)\n               - (c_i*Δz_i/Δt) * (T_f - T_i_n)\n\nWhen is_top = 0.0 exactly, the top-only term vanishes and the interior subtraction is active, exactly reproducing the upstream `:interior` equation; when is_top = 1.0 exactly, the interior subtraction vanishes and the surface-flux term is active, exactly reproducing the upstream `:top` equation. The masking is numerically exact in IEEE-754 (0.0*x = 0.0 for finite x and 1.0*x = x), so neither variant carries any roundoff penalty relative to the upstream non-merged form. The dead parameters in each mode (h_n, dh_dT when is_top=0; F_im1_n, F_im1_np1 when is_top=1) are still declared but multiplied by 0 — they act as null inputs and do not affect the result.\n\nGreek `α_CN` is renamed `alpha_CN` for ASCII parameter naming per the rig convention (matches the heating_cooling_flux.esm sibling's encoding of the same upstream `@constants` α_CN). The other Unicode names (Δz_i, Δt) are preserved verbatim to match sibling ESM files (layer_phase_change_energy.esm uses `Δt`, exponential_grid.esm uses `dz_layer_i`, etc.).\n\nUpstream `L_f = 3.337e5 J/kg` is declared in the `@constants` block of the source `PhaseChangeEnergy` component but does NOT appear in either variant's `eqs` (verified by inspection of the source at :522-:547). It is the latent heat of fusion held for downstream use by sibling components LayerPhaseChangeEnergy (Eq. 4.73) and SnowMeltNoLayers (Eq. 4.71). It is omitted here since the migrated ESM should reflect what the equations actually depend on; the sibling files carry their own L_f literal at the same value (3.337e5 J/kg, Oleson Table 1.4) for the equations that do use it. This mirrors the upstream design where each component re-declares the constants it uses rather than sharing a registry.\n\nObservability classification: this component is functionally observed-only — the single output variable H_excess is given by a single algebraic equation in 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, layer_phase_change_energy.esm, total_phase_change_energy.esm).\n\nTolerance rationale (model-level rel=1e-9): observed-only single-equation algebraic system (no ODE integration); the equation is a sum of products with one division — IEEE-754 ULP precision (~1e-16). 1e-9 leaves comfortable headroom above ULP roundoff and matches the convention used by sibling observed-only files (heating_cooling_flux.esm uses 1e-9 for the analogous Crank-Nicholson HVAC equations, layer_phase_change_energy.esm uses 1e-9 for the analogous Eq. 4.73 single-equation form). Tighter `abs: 1e-12` is asserted on outputs that are exactly zero by construction so the assertion does not silently accept a stuck-at-1e-3 regression.\n\nUpstream tolerance reproduction: upstream `PhaseChangeEnergy Interior` test (test/roof_wall_road_snow_temperatures_test.jl:413) and `PhaseChangeEnergy Top` test (:443) assert `rtol = 1.0e-10`. Re-asserted at the same tolerance via per-assertion `tolerance: { rel: 1.0e-10 }` on the upstream-canonical tests, tighter than the model-level default."
  },
  "tolerance": {
    "rel": 1e-09
  },
  "variables": {
    "T_f": {
      "type": "parameter",
      "units": "K",
      "default": 273.15,
      "description": "Freezing temperature of water (Oleson Table 1.4) — upstream `@constants` T_f; threshold the phase-change energy diagnostic compares the layer temperature against"
    },
    "alpha_CN": {
      "type": "parameter",
      "units": "1",
      "default": 0.5,
      "description": "Crank-Nicholson weight (dimensionless) — upstream `@constants` α_CN; α=0.5 gives the symmetric time-centered average of F_i_n and F_i_np1"
    },
    "c_i": {
      "type": "parameter",
      "units": "J/(m^3*K)",
      "default": 2000000.0,
      "description": "Volumetric heat capacity of layer i (Oleson Eq. 4.59 storage coefficient) — typical soil/snow ~2e6 J/(m^3*K)"
    },
    "Δz_i": {
      "type": "parameter",
      "units": "m",
      "default": 0.1,
      "description": "Layer thickness (Oleson Eq. 4.59) — distance between adjacent interfaces of the thermal column"
    },
    "Δt": {
      "type": "parameter",
      "units": "s",
      "default": 3600.0,
      "description": "Time step (Oleson Eq. 4.59) — implicit-explicit integrator step over which the Crank-Nicholson average is taken"
    },
    "T_i_n": {
      "type": "parameter",
      "units": "K",
      "default": 273.15,
      "description": "Layer temperature at time n (Oleson Eq. 4.59) — input from the urban-canopy temperature solution at the start of the step"
    },
    "F_i_n": {
      "type": "parameter",
      "units": "W/m^2",
      "default": 0.0,
      "description": "Heat flux at the bottom interface of layer i at time n (Oleson Eq. 4.59 lower-boundary Crank-Nicholson contribution); positive = downward into layer i+1"
    },
    "F_i_np1": {
      "type": "parameter",
      "units": "W/m^2",
      "default": 0.0,
      "description": "Heat flux at the bottom interface of layer i at time n+1 (Oleson Eq. 4.59) — paired with F_i_n in the Crank-Nicholson average"
    },
    "F_im1_n": {
      "type": "parameter",
      "units": "W/m^2",
      "default": 0.0,
      "description": "Heat flux at the top interface of layer i (= bottom of layer i-1) at time n (Oleson Eq. 4.59 interior-variant upper-boundary Crank-Nicholson contribution). Active only when is_top = 0.0; multiplied by zero when is_top = 1.0 (top variant uses h_n / dh_dT instead)."
    },
    "F_im1_np1": {
      "type": "parameter",
      "units": "W/m^2",
      "default": 0.0,
      "description": "Heat flux at the top interface of layer i at time n+1 (Oleson Eq. 4.59 interior variant) — paired with F_im1_n in the Crank-Nicholson average. Active only when is_top = 0.0."
    },
    "h_n": {
      "type": "parameter",
      "units": "W/m^2",
      "default": 0.0,
      "description": "Surface heat flux at time n (Oleson Eq. 4.59 top-variant upper-boundary forcing). Active only when is_top = 1.0; multiplied by zero when is_top = 0.0 (interior variant uses F_im1_n / F_im1_np1 instead)."
    },
    "dh_dT": {
      "type": "parameter",
      "units": "W/(m^2*K)",
      "default": 0.0,
      "description": "Linearization slope of the surface heat flux with respect to layer-1 temperature (Oleson Eq. 4.59 top-variant ∂h/∂T term, evaluated at T_f via the (T_f - T_i_n) factor). Active only when is_top = 1.0."
    },
    "is_top": {
      "type": "parameter",
      "units": "1",
      "default": 0.0,
      "description": "Dimensionless switch selecting the upper-boundary variant of Eq. 4.59: 0.0 selects the interior variant (i > snl+1, upper-interface Crank-Nicholson flux), 1.0 selects the top variant (i = snl+1, surface heat flux with linearization). Default 0.0 matches the upstream `layer_type = :interior` default. Other values are not physically meaningful — the merged equation is exact only at the two discrete settings 0.0 and 1.0."
    },
    "H_excess": {
      "type": "unknown",
      "units": "W/m^2",
      "description": "Energy excess (positive) or deficit (negative) available for phase change in layer i (Oleson Eq. 4.59). Computed as net interface-flux divergence minus the storage required to bring the layer to T_f, with the top-vs-interior upper-boundary variant selected by is_top."
    }
  },
  "equations": [
    {
      "lhs": "H_excess",
      "rhs": {
        "op": "-",
        "args": [
          {
            "op": "+",
            "args": [
              {
                "op": "+",
                "args": [
                  {
                    "op": "*",
                    "args": [
                      "alpha_CN",
                      "F_i_n"
                    ]
                  },
                  {
                    "op": "*",
                    "args": [
                      {
                        "op": "-",
                        "args": [
                          1.0,
                          "alpha_CN"
                        ]
                      },
                      "F_i_np1"
                    ]
                  }
                ]
              },
              {
                "op": "*",
                "args": [
                  "is_top",
                  {
                    "op": "+",
                    "args": [
                      "h_n",
                      {
                        "op": "*",
                        "args": [
                          "dh_dT",
                          {
                            "op": "-",
                            "args": [
                              "T_f",
                              "T_i_n"
                            ]
                          }
                        ]
                      }
                    ]
                  }
                ]
              }
            ]
          },
          {
            "op": "+",
            "args": [
              {
                "op": "*",
                "args": [
                  {
                    "op": "-",
                    "args": [
                      1.0,
                      "is_top"
                    ]
                  },
                  {
                    "op": "+",
                    "args": [
                      {
                        "op": "*",
                        "args": [
                          "alpha_CN",
                          "F_im1_n"
                        ]
                      },
                      {
                        "op": "*",
                        "args": [
                          {
                            "op": "-",
                            "args": [
                              1.0,
                              "alpha_CN"
                            ]
                          },
                          "F_im1_np1"
                        ]
                      }
                    ]
                  }
                ]
              },
              {
                "op": "*",
                "args": [
                  {
                    "op": "/",
                    "args": [
                      {
                        "op": "*",
                        "args": [
                          "c_i",
                          "Δz_i"
                        ]
                      },
                      "Δt"
                    ]
                  },
                  {
                    "op": "-",
                    "args": [
                      "T_f",
                      "T_i_n"
                    ]
                  }
                ]
              }
            ]
          }
        ]
      }
    }
  ],
  "tests": [
    {
      "id": "upstream_interior_canonical",
      "description": "Reproduces upstream `PhaseChangeEnergy Interior` test (UrbanCanopy.jl test/roof_wall_road_snow_temperatures_test.jl:413). is_top = 0.0 selects the interior variant. Inputs: c_i = 2e6 J/(m^3*K), Δz_i = 0.1 m, Δt = 3600 s, T_i_n = 274 K (above T_f = 273.15 K → negative storage term as the layer cools toward freezing); F_i_n = -10, F_im1_n = -5, F_i_np1 = -12, F_im1_np1 = -6 W/m^2 (downward fluxes; flux divergence -5 → -6 W/m^2 at time n/n+1). Expected H_excess = α*(F_i_n - F_im1_n) + (1-α)*(F_i_np1 - F_im1_np1) - (c_i*Δz_i/Δt)*(T_f - T_i_n) = 0.5*(-5) + 0.5*(-6) - (2e6*0.1/3600)*(-0.85) = -2.5 - 3.0 + 47.222... = 41.7222 W/m^2 (positive: storage gain dominates, layer can release energy to phase change). rtol matches upstream `rtol = 1.0e-10`.",
      "parameter_overrides": {
        "is_top": 0.0,
        "c_i": 2000000.0,
        "Δz_i": 0.1,
        "Δt": 3600.0,
        "T_i_n": 274.0,
        "F_i_n": -10.0,
        "F_i_np1": -12.0,
        "F_im1_n": -5.0,
        "F_im1_np1": -6.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": 41.72222222222222,
          "tolerance": {
            "rel": 1e-10
          }
        }
      ]
    },
    {
      "id": "upstream_top_canonical",
      "description": "Reproduces upstream `PhaseChangeEnergy Top` test (UrbanCanopy.jl test/roof_wall_road_snow_temperatures_test.jl:443). is_top = 1.0 selects the top variant. Inputs: c_i = 2e6, Δz_i = 0.01 m (thin top layer), Δt = 3600 s, T_i_n = 272 K (below T_f → positive storage term as the layer warms toward freezing); h_n = 100 W/m^2 (incoming surface flux), dh_dT = -15 W/(m^2*K) (linearization slope, negative: more flux when the layer is colder); F_i_n = -5, F_i_np1 = -6 W/m^2. Expected H_excess = h_n + dh_dT*(T_f - T_i_n) + α*F_i_n + (1-α)*F_i_np1 - (c_i*Δz_i/Δt)*(T_f - T_i_n) = 100 + (-15)*1.15 + 0.5*(-5) + 0.5*(-6) - (2e6*0.01/3600)*1.15 = 100 - 17.25 - 2.5 - 3.0 - 6.3888... = 70.8611 W/m^2 (positive: surface forcing dominates). rtol matches upstream `rtol = 1.0e-10`.",
      "parameter_overrides": {
        "is_top": 1.0,
        "c_i": 2000000.0,
        "Δz_i": 0.01,
        "Δt": 3600.0,
        "T_i_n": 272.0,
        "h_n": 100.0,
        "dh_dT": -15.0,
        "F_i_n": -5.0,
        "F_i_np1": -6.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": 70.86111111111111,
          "tolerance": {
            "rel": 1e-10
          }
        }
      ]
    },
    {
      "id": "interior_at_freezing_no_storage",
      "description": "Interior variant boundary case: T_i_n = T_f = 273.15 K exactly so the storage term (c_i*Δz_i/Δt)*(T_f - T_i_n) is exactly zero by construction. Only the Crank-Nicholson flux divergence contributes: H_excess = α*(F_i_n - F_im1_n) + (1-α)*(F_i_np1 - F_im1_np1) = 0.5*(20 - 5) + 0.5*(15 - 10) = 7.5 + 2.5 = 10.0 W/m^2 (positive divergence = energy accumulating in this layer because more flows in from below than flows up to layer i-1). Isolates the flux-divergence term from the storage term.",
      "parameter_overrides": {
        "is_top": 0.0,
        "c_i": 2000000.0,
        "Δz_i": 0.1,
        "Δt": 3600.0,
        "T_i_n": 273.15,
        "F_i_n": 20.0,
        "F_i_np1": 15.0,
        "F_im1_n": 5.0,
        "F_im1_np1": 10.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": 10.0
        }
      ]
    },
    {
      "id": "interior_balanced_fluxes_zero_excess",
      "description": "Interior variant equilibrium case: T_i_n = T_f exactly AND F_i_n = F_im1_n = F_i_np1 = F_im1_np1 = 8 W/m^2 (same downward flux through both interfaces, both time levels). Both Crank-Nicholson differences vanish (8 - 8 = 0) and the storage term vanishes (T_f - T_i_n = 0). H_excess = 0 — no surplus or deficit, the layer is in thermal equilibrium. Exercises the all-zero result of the merged equation when both contributions independently cancel.",
      "parameter_overrides": {
        "is_top": 0.0,
        "c_i": 2000000.0,
        "Δz_i": 0.1,
        "Δt": 3600.0,
        "T_i_n": 273.15,
        "F_i_n": 8.0,
        "F_i_np1": 8.0,
        "F_im1_n": 8.0,
        "F_im1_np1": 8.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": 0.0,
          "tolerance": {
            "abs": 1e-12
          }
        }
      ]
    },
    {
      "id": "interior_cold_layer_positive_storage",
      "description": "Interior variant with T_i_n = 263.15 K (10 K below T_f); positive storage term (c_i*Δz_i/Δt)*(T_f - T_i_n) = (2e6*0.1/3600)*10 = 555.555 W/m^2 would be subtracted, meaning the layer needs to ABSORB 555.555 W/m^2 just to warm to T_f. With balanced fluxes (F_i_n = F_im1_n = F_i_np1 = F_im1_np1 = 0), the Crank-Nicholson terms are zero and H_excess = -555.555 W/m^2 — large energy deficit. Anchors the sign convention: H_excess < 0 means the layer cannot reach T_f within Δt without additional energy.",
      "parameter_overrides": {
        "is_top": 0.0,
        "c_i": 2000000.0,
        "Δz_i": 0.1,
        "Δt": 3600.0,
        "T_i_n": 263.15,
        "F_i_n": 0.0,
        "F_i_np1": 0.0,
        "F_im1_n": 0.0,
        "F_im1_np1": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": -555.5555555555555,
          "tolerance": {
            "rel": 1e-10
          }
        }
      ]
    },
    {
      "id": "interior_warm_layer_negative_storage",
      "description": "Mirror of `interior_cold_layer_positive_storage`: T_i_n = 283.15 K (10 K above T_f); the storage term flips sign. (c_i*Δz_i/Δt)*(T_f - T_i_n) = (2e6*0.1/3600)*(-10) = -555.555 W/m^2 — subtracted to give +555.555 W/m^2 surplus. With all fluxes zero, H_excess = +555.555 W/m^2 — large positive excess available to melt ice if any is present. Confirms the storage term is antisymmetric in (T_f - T_i_n) around T_f.",
      "parameter_overrides": {
        "is_top": 0.0,
        "c_i": 2000000.0,
        "Δz_i": 0.1,
        "Δt": 3600.0,
        "T_i_n": 283.15,
        "F_i_n": 0.0,
        "F_i_np1": 0.0,
        "F_im1_n": 0.0,
        "F_im1_np1": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": 555.5555555555555,
          "tolerance": {
            "rel": 1e-10
          }
        }
      ]
    },
    {
      "id": "top_no_surface_flux_reduces_to_bottom_only",
      "description": "Top variant boundary case: h_n = 0 and dh_dT = 0 → the surface-flux term vanishes; T_i_n = T_f → storage term vanishes. With F_i_n = F_i_np1 = 50 W/m^2 (uniform downward bottom flux), only the Crank-Nicholson bottom-flux term contributes: H_excess = α*F_i_n + (1-α)*F_i_np1 + 0 - 0 = 0.5*50 + 0.5*50 = 50.0 W/m^2. Exercises the degenerate top variant with zero surface forcing — it should produce the bare bottom-flux contribution without contamination from the (multiplied-by-zero) F_im1 terms.",
      "parameter_overrides": {
        "is_top": 1.0,
        "c_i": 2000000.0,
        "Δz_i": 0.01,
        "Δt": 3600.0,
        "T_i_n": 273.15,
        "h_n": 0.0,
        "dh_dT": 0.0,
        "F_i_n": 50.0,
        "F_i_np1": 50.0,
        "F_im1_n": 999.0,
        "F_im1_np1": -999.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": 50.0
        }
      ]
    },
    {
      "id": "top_dead_F_im1_inputs_ignored",
      "description": "Top variant masking-check: is_top = 1.0 so the (1 - is_top) factor on the interior-variant F_im1 terms is exactly zero — even huge F_im1_n = F_im1_np1 = 1e8 W/m^2 (unphysical magnitude) must not affect H_excess. With T_i_n = T_f, h_n = 100, dh_dT = 0, F_i_n = F_i_np1 = 0: H_excess should equal h_n = 100 W/m^2. If the masking is wrong (e.g. interior subtraction not gated), H_excess would be off by 0.5e8 + 0.5e8 = 1e8 — a 1e6× error easy to detect. Confirms the is_top switch is exact at 1.0.",
      "parameter_overrides": {
        "is_top": 1.0,
        "c_i": 2000000.0,
        "Δz_i": 0.01,
        "Δt": 3600.0,
        "T_i_n": 273.15,
        "h_n": 100.0,
        "dh_dT": 0.0,
        "F_i_n": 0.0,
        "F_i_np1": 0.0,
        "F_im1_n": 100000000.0,
        "F_im1_np1": 100000000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": 100.0
        }
      ]
    },
    {
      "id": "interior_dead_h_dh_inputs_ignored",
      "description": "Mirror of `top_dead_F_im1_inputs_ignored`: is_top = 0.0 so the is_top factor on the top-variant (h_n + dh_dT*(T_f - T_i_n)) term is exactly zero — even huge h_n = 1e8 W/m^2 and dh_dT = 1e6 W/(m^2*K) must not affect H_excess. With T_i_n = T_f, all bottom and upper-interface fluxes balanced at 0: H_excess = 0. Confirms the is_top switch is exact at 0.0 and the upstream interior variant is recovered when the top-variant inputs are noise.",
      "parameter_overrides": {
        "is_top": 0.0,
        "c_i": 2000000.0,
        "Δz_i": 0.1,
        "Δt": 3600.0,
        "T_i_n": 273.15,
        "h_n": 100000000.0,
        "dh_dT": 1000000.0,
        "F_i_n": 0.0,
        "F_i_np1": 0.0,
        "F_im1_n": 0.0,
        "F_im1_np1": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": 0.0,
          "tolerance": {
            "abs": 1e-12
          }
        }
      ]
    },
    {
      "id": "top_dh_dT_linearization_amplifies",
      "description": "Top variant: probes the (T_f - T_i_n) linearization slope of the surface flux. T_i_n = 268.15 K (5 K below T_f), h_n = 0, dh_dT = -20 W/(m^2*K) — when the layer is cold, the surface flux is enhanced relative to its T_f linearization point. Surface-flux term = h_n + dh_dT*(T_f - T_i_n) = 0 + (-20)*5 = -100 W/m^2. Storage term = (2e6*0.01/3600)*5 = 27.778 W/m^2 — subtracted. With F_i_n = F_i_np1 = 0: H_excess = 0 + (-100) - 0 - 27.778 = -127.778 W/m^2. Exercises the dh_dT contribution sign and magnitude.",
      "parameter_overrides": {
        "is_top": 1.0,
        "c_i": 2000000.0,
        "Δz_i": 0.01,
        "Δt": 3600.0,
        "T_i_n": 268.15,
        "h_n": 0.0,
        "dh_dT": -20.0,
        "F_i_n": 0.0,
        "F_i_np1": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": -127.77777777777777,
          "tolerance": {
            "rel": 1e-10
          }
        }
      ]
    },
    {
      "id": "interior_asymmetric_crank_nicholson",
      "description": "Interior variant: stress-tests the Crank-Nicholson averaging with asymmetric inputs at times n vs n+1. T_i_n = T_f, F_i_n = 100, F_im1_n = 0, F_i_np1 = -100, F_im1_np1 = 0 (sign-flipped bottom flux during the step; upper interface inert). H_excess = α*(F_i_n - F_im1_n) + (1-α)*(F_i_np1 - F_im1_np1) - 0 = 0.5*100 + 0.5*(-100) = 0 — the Crank-Nicholson combiner cancels the sign flip. Verifies that the time-averaging is symmetric in (F_*_n, F_*_np1).",
      "parameter_overrides": {
        "is_top": 0.0,
        "c_i": 2000000.0,
        "Δz_i": 0.1,
        "Δt": 3600.0,
        "T_i_n": 273.15,
        "F_i_n": 100.0,
        "F_i_np1": -100.0,
        "F_im1_n": 0.0,
        "F_im1_np1": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": 0.0,
          "tolerance": {
            "abs": 1e-12
          }
        }
      ]
    },
    {
      "id": "interior_thin_layer_short_step_scaling",
      "description": "Interior variant: probes the storage-coefficient scaling (c_i*Δz_i/Δt). Use a thin layer (Δz_i = 0.001 m, ~1 mm) and a short step (Δt = 60 s, ~1 min) with a 1 K temperature deviation (T_i_n = 274.15, T_f - T_i_n = -1 K). Storage = (2e6*0.001/60)*(-1) = -33.333... W/m^2 subtracted → +33.333 W/m^2 surplus. With all fluxes zero: H_excess = +33.333 W/m^2. Confirms the storage term scales correctly with the thin-layer, short-step thermal-inertia limit (large c_i*Δz_i/Δt with small T deviation).",
      "parameter_overrides": {
        "is_top": 0.0,
        "c_i": 2000000.0,
        "Δz_i": 0.001,
        "Δt": 60.0,
        "T_i_n": 274.15,
        "F_i_n": 0.0,
        "F_i_np1": 0.0,
        "F_im1_n": 0.0,
        "F_im1_np1": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "H_excess",
          "time": 0.0,
          "expected": 33.33333333333333,
          "tolerance": {
            "rel": 1e-10
          }
        }
      ]
    }
  ],
  "analyses": [
    {
      "id": "H_excess_vs_T_i_n_interior",
      "description": "Interior variant: sweep T_i_n across the freezing point T_f = 273.15 K with all fluxes zero. With F_i_n = F_im1_n = F_i_np1 = F_im1_np1 = 0, only the storage term contributes: H_excess = -(c_i*Δz_i/Δt)*(T_f - T_i_n). At T_i_n = T_f the curve crosses zero; for T_i_n > T_f the layer has surplus energy (warm layer can melt ice if any); for T_i_n < T_f the layer has deficit (cold layer cannot reach T_f within Δt without additional energy). Linear in T_i_n with slope +(c_i*Δz_i/Δt) = 55.555 W/(m^2*K) at the default parameters. Surfaces the Eq. 4.59 storage-term sign convention directly.",
      "initial_state": {},
      "parameters": {
        "is_top": 0.0,
        "c_i": 2000000.0,
        "Δz_i": 0.1,
        "Δt": 3600.0,
        "F_i_n": 0.0,
        "F_i_np1": 0.0,
        "F_im1_n": 0.0,
        "F_im1_np1": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "parameter_sweep": {
        "type": "cartesian",
        "dimensions": [
          {
            "parameter": "T_i_n",
            "range": {
              "start": 263.15,
              "stop": 283.15,
              "count": 41,
              "scale": "linear"
            }
          }
        ]
      },
      "plots": [
        {
          "id": "H_excess_curve_interior",
          "type": "line",
          "description": "H_excess vs T_i_n (interior, zero fluxes): straight line through (T_f, 0) with slope c_i*Δz_i/Δt — energy deficit for cold layers, energy surplus for warm layers.",
          "x": {
            "variable": "T_i_n",
            "label": "T_i_n (K)"
          },
          "y": {
            "variable": "H_excess",
            "label": "H_excess (W/m^2)"
          }
        }
      ]
    },
    {
      "id": "H_excess_vs_T_i_n_top",
      "description": "Top variant: same temperature sweep as the interior example but with is_top = 1.0 and a representative surface forcing h_n = 100 W/m^2 with dh_dT = -10 W/(m^2*K). The top variant adds h_n + dh_dT*(T_f - T_i_n) to the layer's energy balance — when the layer is cold (T_i_n < T_f), the linearization slope amplifies the surface flux (negative dh_dT × positive (T_f - T_i_n) > 0 → enhanced flux into the layer when cold); when the layer is warm, the flux is reduced. Compared to the interior example, the curve is shifted upward by the surface-forcing contribution and tilted by the dh_dT slope.",
      "initial_state": {},
      "parameters": {
        "is_top": 1.0,
        "c_i": 2000000.0,
        "Δz_i": 0.01,
        "Δt": 3600.0,
        "h_n": 100.0,
        "dh_dT": -10.0,
        "F_i_n": 0.0,
        "F_i_np1": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "parameter_sweep": {
        "type": "cartesian",
        "dimensions": [
          {
            "parameter": "T_i_n",
            "range": {
              "start": 263.15,
              "stop": 283.15,
              "count": 41,
              "scale": "linear"
            }
          }
        ]
      },
      "plots": [
        {
          "id": "H_excess_curve_top",
          "type": "line",
          "description": "H_excess vs T_i_n (top variant, surface forcing h_n = 100, dh_dT = -10): tilted line offset upward by h_n with slope c_i*Δz_i/Δt - dh_dT relative to the interior example.",
          "x": {
            "variable": "T_i_n",
            "label": "T_i_n (K)"
          },
          "y": {
            "variable": "H_excess",
            "label": "H_excess (W/m^2)"
          }
        }
      ]
    },
    {
      "id": "two_d_sweep_T_F_interior",
      "description": "Interior variant: 2-D sweep over (T_i_n, F_i_n) at fixed F_im1_n = F_i_np1 = F_im1_np1 = 0 surfacing the joint contribution of the storage term (linear in T_i_n) and the Crank-Nicholson lower-interface flux (linear in F_i_n). Iso-H_excess contours are oblique lines in (T_i_n, F_i_n) space with slope determined by the storage coefficient and α_CN. Exercises the bivariate dependence of Eq. 4.59 with the upper-interface fluxes set to zero (limit of a column with no flux at the layer-(i-1)/layer-i interface).",
      "initial_state": {},
      "parameters": {
        "is_top": 0.0,
        "c_i": 2000000.0,
        "Δz_i": 0.1,
        "Δt": 3600.0,
        "F_i_np1": 0.0,
        "F_im1_n": 0.0,
        "F_im1_np1": 0.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "parameter_sweep": {
        "type": "cartesian",
        "dimensions": [
          {
            "parameter": "T_i_n",
            "range": {
              "start": 263.15,
              "stop": 283.15,
              "count": 21,
              "scale": "linear"
            }
          },
          {
            "parameter": "F_i_n",
            "range": {
              "start": -100.0,
              "stop": 100.0,
              "count": 21,
              "scale": "linear"
            }
          }
        ]
      },
      "plots": [
        {
          "id": "H_excess_T_F_heatmap",
          "type": "heatmap",
          "description": "H_excess in (T_i_n, F_i_n) space (interior variant). Planar surface: storage slope along T_i_n and α_CN-weighted Crank-Nicholson slope along F_i_n. Iso-H_excess lines have constant slope, with the zero contour passing through (T_f, 0).",
          "x": {
            "variable": "T_i_n",
            "label": "T_i_n (K)"
          },
          "y": {
            "variable": "F_i_n",
            "label": "F_i_n (W/m^2)"
          },
          "value": {
            "variable": "H_excess",
            "reduce": "final"
          }
        }
      ]
    }
  ]
}