FreezingPointDepression

{
  "reference": {
    "notes": "source_ref: UrbanCanopy.jl@9200edd4f15609f20339d3a7d104de8064544310:src/roof_wall_road_snow_temperatures.jl:949\nUpstream `FreezingPointDepression` is an MTK `@component` whose six inputs (T_i, θ_sat, Δz, ψ_sat, B_i, ρ_liq) are `@parameters` and whose single output (w_liq_max) is a time-indexed `@variable` with one algebraic equation (Eq. 4.58 with an `ifelse` branch on T_i < T_f). Identical equation structure preserved here.\n\nThe upstream component declares four `@constants`: T_f (273.15 K), L_f (3.337e5 J/kg), g (9.80616 m/s^2), and one_ref (1.0, dimensionless reference for non-dimensionalization of the bracket prior to the fractional exponent). They are encoded here as parameters with the upstream values as `default`. one_ref is kept in the equation to maintain 1:1 fidelity with the upstream MTK unit-balancing pattern (the same convention used by snow_melt_no_layers.esm for its `zero_kgm2`/`zero_m` literals): the fractional exponent in Eq. 4.58 is well-defined only on a dimensionless argument, and one_ref makes the dimensional cancellation explicit to MTK's unit checker upstream.\n\nUnicode parameter names (Δz, θ_sat, ψ_sat, ρ_liq) are preserved from upstream — peer ESM components in this rig (snow_melt_no_layers.esm uses Δt, snow_density.esm uses ρ_*) already use unicode in identifiers and the schema's identifier pattern allows it.\n\nThe upstream comment maps the equation to CLM Tech Note Eq. 4.58 and notes a unit subtlety: the paper uses ψ_sat in mm but the upstream component takes ψ_sat in SI (m); the 10^3 mm→m conversion factor is absorbed (so no factor-of-1000 appears in the equation). ψ_sat is conventionally negative (soil matric potential); the equation negates it to keep the base of the fractional exponent positive. These are preserved in the equation as-is.\n\nObservability classification: this component is functionally observed-only — the single output variable is given by a single algebraic equation in terms of the parameters (the `(t)` time index has no role since there is no derivative). 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 (snow_melt_no_layers.esm, total_phase_change_energy.esm, ideal_gas.esm).\n\nTolerance rationale (model-level rel=1e-9): observed-only single-equation algebraic system (no ODE integration); the Eq. 4.58 evaluation is a chain of multiplies/divides plus one fractional-power call. Roundoff is at the ULP level (~1e-16) for the multiplicative parts and ~1e-15 for `pow`. 1e-9 leaves comfortable headroom and matches sibling temps files (total_phase_change_energy.esm uses 1e-9). Tighter `abs: 1e-12` is asserted on the `ifelse`-else branch outputs (T_i > T_f cases) where the value reduces to the exact product Δz·θ_sat·ρ_liq with no fractional-power term — the assertion would not silently accept a stuck-at-1e-3 regression. Reference values one-shot generated from the canonical formula (identical to upstream MTK eq.) to IEEE-754 double precision.\n\nThe `T_i = T_f` exact-boundary cannot be tested through the canonical Python runner: the runner's sympy.Piecewise eagerly evaluates both branches symbolically, and at T_i = T_f the cold branch's bracket is 0 raised to the negative power -1/B_i — undefined. The upstream Julia MTK runtime ifelse only evaluates the selected branch and is unaffected. To still anchor the boundary convention (the upstream uses strict `<`, so equality belongs to the saturated `else` branch), the test pair `epsilon_below_freezing` (T_i = 273.149, cold branch, 191.4 kg/m^2) and `just_above_freezing` (T_i = 273.151, else branch, exactly 200 kg/m^2) brackets the boundary at the 1-mK scale — tight enough to detect any swap of `<` and `≤`."
  },
  "tolerance": {
    "rel": 1e-09
  },
  "variables": {
    "T_i": {
      "type": "parameter",
      "units": "K",
      "default": 280.0,
      "description": "Layer temperature"
    },
    "θ_sat": {
      "type": "parameter",
      "units": "1",
      "default": 0.4,
      "description": "Porosity (dimensionless)"
    },
    "Δz": {
      "type": "parameter",
      "units": "m",
      "default": 0.5,
      "description": "Layer thickness"
    },
    "ψ_sat": {
      "type": "parameter",
      "units": "m",
      "default": -0.1,
      "description": "Saturated matric potential (negative by convention; SI m, not mm)"
    },
    "B_i": {
      "type": "parameter",
      "units": "1",
      "default": 5.0,
      "description": "Clapp and Hornberger exponent (dimensionless)"
    },
    "ρ_liq": {
      "type": "parameter",
      "units": "kg/m^3",
      "default": 1000.0,
      "description": "Density of liquid water"
    },
    "T_f": {
      "type": "parameter",
      "units": "K",
      "default": 273.15,
      "description": "Freezing temperature of water (Oleson et al. 2010 Table 1.4 — 273.15 K)"
    },
    "L_f": {
      "type": "parameter",
      "units": "J/kg",
      "default": 333700.0,
      "description": "Latent heat of fusion of water (Oleson et al. 2010 Table 1.4 — 3.337e5 J/kg)"
    },
    "g": {
      "type": "parameter",
      "units": "m/s^2",
      "default": 9.80616,
      "description": "Gravitational acceleration (Oleson et al. 2010 Table 1.4 — 9.80616 m/s^2)"
    },
    "one_ref": {
      "type": "parameter",
      "units": "1",
      "default": 1.0,
      "description": "Dimensionless reference for non-dimensionalization of the bracket in Eq. 4.58 prior to the fractional exponent — preserved 1:1 from upstream MTK unit-balancing pattern"
    },
    "w_liq_max": {
      "type": "unknown",
      "units": "kg/m^2",
      "description": "Maximum liquid water content when T_i < T_f (Oleson Eq. 4.58); above T_f equals fully-saturated Δz·θ_sat·ρ_liq"
    }
  },
  "equations": [
    {
      "lhs": "w_liq_max",
      "rhs": {
        "op": "ifelse",
        "args": [
          {
            "op": "<",
            "args": [
              "T_i",
              "T_f"
            ]
          },
          {
            "op": "*",
            "args": [
              {
                "op": "*",
                "args": [
                  {
                    "op": "*",
                    "args": [
                      "Δz",
                      "θ_sat"
                    ]
                  },
                  "ρ_liq"
                ]
              },
              {
                "op": "^",
                "args": [
                  {
                    "op": "*",
                    "args": [
                      {
                        "op": "/",
                        "args": [
                          {
                            "op": "*",
                            "args": [
                              "L_f",
                              {
                                "op": "-",
                                "args": [
                                  "T_f",
                                  "T_i"
                                ]
                              }
                            ]
                          },
                          {
                            "op": "*",
                            "args": [
                              {
                                "op": "*",
                                "args": [
                                  "g",
                                  "T_i"
                                ]
                              },
                              {
                                "op": "-",
                                "args": [
                                  "ψ_sat"
                                ]
                              }
                            ]
                          }
                        ]
                      },
                      "one_ref"
                    ]
                  },
                  {
                    "op": "/",
                    "args": [
                      -1,
                      "B_i"
                    ]
                  }
                ]
              }
            ]
          },
          {
            "op": "*",
            "args": [
              {
                "op": "*",
                "args": [
                  "Δz",
                  "θ_sat"
                ]
              },
              "ρ_liq"
            ]
          }
        ]
      }
    }
  ],
  "tests": [
    {
      "id": "above_freezing_saturated",
      "description": "Reproduces upstream `FreezingPointDepression - Above Freezing` test (UrbanCanopy.jl test/roof_wall_road_snow_temperatures_test.jl:1052). T_i = 280 K (above T_f), θ_sat = 0.4, Δz = 0.5 m, ψ_sat = -0.1 m, B_i = 5.0, ρ_liq = 1000 kg/m^3 → w_liq_max = Δz·θ_sat·ρ_liq = 200 kg/m^2 exactly (else branch of the ifelse, no fractional-power evaluation).",
      "parameter_overrides": {
        "T_i": 280.0,
        "θ_sat": 0.4,
        "Δz": 0.5,
        "ψ_sat": -0.1,
        "B_i": 5.0,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "w_liq_max",
          "time": 0.0,
          "expected": 200.0,
          "tolerance": {
            "abs": 1e-12
          }
        }
      ]
    },
    {
      "id": "just_above_freezing",
      "description": "T_i = 273.151 K (1 mK above T_f) — the just-above-freezing side of the boundary. Falls into the `else` branch (saturated max = 200 kg/m^2). Paired with `epsilon_below_freezing` (1 mK on the other side, 191.4 kg/m^2) to anchor the discontinuity at T_f imposed by the upstream `ifelse(T_i < T_f, …)`. (T_i = T_f exactly cannot be tested through the canonical Python runner — the runner's sympy.Piecewise eagerly evaluates both branches symbolically, and the cold branch's `(0)^(-1/B_i)` is undefined; the upstream Julia MTK ifelse only evaluates the selected branch at runtime, so the upstream is unaffected. The 1-mK pair brackets the boundary tightly enough.)",
      "parameter_overrides": {
        "T_i": 273.151,
        "θ_sat": 0.4,
        "Δz": 0.5,
        "ψ_sat": -0.1,
        "B_i": 5.0,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "w_liq_max",
          "time": 0.0,
          "expected": 200.0,
          "tolerance": {
            "abs": 1e-12
          }
        }
      ]
    },
    {
      "id": "below_freezing_5K",
      "description": "Reproduces upstream `FreezingPointDepression - Below Freezing` test (UrbanCanopy.jl test/roof_wall_road_snow_temperatures_test.jl:1073). T_i = 268.15 K (5 K below T_f), θ_sat = 0.4, Δz = 0.5 m, ψ_sat = -0.1 m, B_i = 5.0, ρ_liq = 1000 kg/m^3. The upstream test only checks 0 < w_liq_max < 200 (max_possible); we anchor the exact value computed from the canonical Eq. 4.58 formula. Verifies the fractional-power branch of the ifelse.",
      "parameter_overrides": {
        "T_i": 268.15,
        "θ_sat": 0.4,
        "Δz": 0.5,
        "ψ_sat": -0.1,
        "B_i": 5.0,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "w_liq_max",
          "time": 0.0,
          "expected": 34.71683367120432
        }
      ]
    },
    {
      "id": "slightly_below_freezing_1K",
      "description": "T_i = 272.15 K (1 K below T_f) — the warm-end of the supercooled-water regime. Anchors the upstream `Colder Means Less Liquid` first-anchor (UrbanCanopy.jl test:1118): w_liq_max at T_i = 272.15 must exceed the value at T_i = 263.15 (paired with `colder_means_less_liquid_anchor` below). Reference value computed from Eq. 4.58.",
      "parameter_overrides": {
        "T_i": 272.15,
        "θ_sat": 0.4,
        "Δz": 0.5,
        "ψ_sat": -0.1,
        "B_i": 5.0,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "w_liq_max",
          "time": 0.0,
          "expected": 48.04190466488627
        }
      ]
    },
    {
      "id": "colder_means_less_liquid_anchor",
      "description": "T_i = 263.15 K (10 K below T_f) — the cold-end anchor of the upstream `Colder Means Less Liquid` test (UrbanCanopy.jl test:1123). Paired with `slightly_below_freezing_1K`: 30.109 < 48.042 confirms the monotone-decreasing relation between w_liq_max and (T_f - T_i) that Niu and Yang (2006) require.",
      "parameter_overrides": {
        "T_i": 263.15,
        "θ_sat": 0.4,
        "Δz": 0.5,
        "ψ_sat": -0.1,
        "B_i": 5.0,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "w_liq_max",
          "time": 0.0,
          "expected": 30.109200535801378
        }
      ]
    },
    {
      "id": "deep_freeze",
      "description": "T_i = 253.15 K (20 K below T_f) — extreme winter regime. The fractional power compresses w_liq_max further to ~26 kg/m^2. Stress-tests the formula in the regime where (T_f - T_i) is large and the bracket inside the (-1/B_i) power is at its largest, ensuring no overflow / loss of precision in the deep-cold corner.",
      "parameter_overrides": {
        "T_i": 253.15,
        "θ_sat": 0.4,
        "Δz": 0.5,
        "ψ_sat": -0.1,
        "B_i": 5.0,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "w_liq_max",
          "time": 0.0,
          "expected": 26.00926841246697
        }
      ]
    },
    {
      "id": "epsilon_below_freezing",
      "description": "T_i = 273.149 K — 1 mK below T_f. (T_f - T_i) is at the small-ΔT scale; w_liq_max is just barely below the saturated max (191.4 vs 200). Verifies the smooth limit of Eq. 4.58 as T_i → T_f^- and that the `<` branch evaluates without spurious cancellation in the small-ΔT regime.",
      "parameter_overrides": {
        "T_i": 273.149,
        "θ_sat": 0.4,
        "Δz": 0.5,
        "ψ_sat": -0.1,
        "B_i": 5.0,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "w_liq_max",
          "time": 0.0,
          "expected": 191.39847463987002
        }
      ]
    },
    {
      "id": "clay_soil_high_B",
      "description": "Clay-rich soil regime: B_i = 10 (high Clapp-Hornberger exponent), θ_sat = 0.5, Δz = 1 m, ψ_sat = -0.3 m at T_i = 263.15 K. The high B_i drives the (-1/B_i) exponent toward 0, so the fractional-power factor approaches 1 and w_liq_max approaches the full saturation Δz·θ_sat·ρ_liq = 500. Result 216.5 < 500 anchors the high-B_i regime.",
      "parameter_overrides": {
        "T_i": 263.15,
        "θ_sat": 0.5,
        "Δz": 1.0,
        "ψ_sat": -0.3,
        "B_i": 10.0,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "w_liq_max",
          "time": 0.0,
          "expected": 216.5293359874589
        }
      ]
    },
    {
      "id": "sandy_soil_low_B_thin",
      "description": "Sandy soil regime: B_i = 3 (low Clapp-Hornberger exponent), θ_sat = 0.3, Δz = 0.05 m (thin), ψ_sat = -0.05 m at T_i = 270 K (3.15 K below T_f). The low B_i drives the (-1/B_i) exponent toward larger negative values, strongly compressing the fractional-power factor; combined with low porosity and thin layer, w_liq_max ~0.75 kg/m^2 is small. Anchors the low-B_i / thin-layer corner.",
      "parameter_overrides": {
        "T_i": 270.0,
        "θ_sat": 0.3,
        "Δz": 0.05,
        "ψ_sat": -0.05,
        "B_i": 3.0,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "assertions": [
        {
          "variable": "w_liq_max",
          "time": 0.0,
          "expected": 0.7518766309616047
        }
      ]
    }
  ],
  "analyses": [
    {
      "id": "w_liq_max_vs_T_i",
      "description": "Sweep T_i across the full supercooled-soil-water regime — from 250 K (deep freeze) through 273.15 (T_f boundary) and into 280 K (above freezing). Reproduces the canonical Niu-Yang (2006) supercooled-water curve: w_liq_max rises monotonically with T_i in the cold branch, hits a discontinuity at T_f (the ifelse boundary), and plateaus at Δz·θ_sat·ρ_liq above freezing.",
      "initial_state": {},
      "parameters": {
        "θ_sat": 0.4,
        "Δz": 0.5,
        "ψ_sat": -0.1,
        "B_i": 5.0,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "parameter_sweep": {
        "type": "cartesian",
        "dimensions": [
          {
            "parameter": "T_i",
            "range": {
              "start": 250.0,
              "stop": 280.0,
              "count": 31,
              "scale": "linear"
            }
          }
        ]
      },
      "plots": [
        {
          "id": "w_liq_max_curve",
          "type": "line",
          "description": "w_liq_max vs T_i: smooth rise in the supercooled-water branch (T_i < T_f), discontinuity at T_f = 273.15 K (ifelse boundary), plateau at 200 kg/m^2 above T_f.",
          "x": {
            "variable": "T_i",
            "label": "T_i (K)"
          },
          "y": {
            "variable": "w_liq_max",
            "label": "w_liq,max (kg/m^2)"
          }
        }
      ]
    },
    {
      "id": "w_liq_max_vs_B_i",
      "description": "Sweep the Clapp-Hornberger exponent B_i across realistic soil-texture range (3 sand → 12 clay) at fixed deep-cold T_i = 263.15 K. Surfaces the role of B_i in compressing w_liq_max via the (-1/B_i) fractional power: low B_i → large compression, high B_i → exponent → 0, w_liq_max → saturated max.",
      "initial_state": {},
      "parameters": {
        "T_i": 263.15,
        "θ_sat": 0.4,
        "Δz": 0.5,
        "ψ_sat": -0.1,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "parameter_sweep": {
        "type": "cartesian",
        "dimensions": [
          {
            "parameter": "B_i",
            "range": {
              "start": 3.0,
              "stop": 12.0,
              "count": 19,
              "scale": "linear"
            }
          }
        ]
      },
      "plots": [
        {
          "id": "w_liq_max_vs_B_curve",
          "type": "line",
          "description": "w_liq_max vs B_i at T_i = 263.15 K: monotone rise toward the saturated max as B_i increases — the (-1/B_i) exponent → 0 limit.",
          "x": {
            "variable": "B_i",
            "label": "B_i (dimensionless)"
          },
          "y": {
            "variable": "w_liq_max",
            "label": "w_liq,max (kg/m^2)"
          }
        }
      ]
    },
    {
      "id": "two_d_sweep_T_B",
      "description": "Joint sweep over T_i (cold branch only, 250-273 K) and B_i (sand-to-clay, 3-12) at fixed θ_sat = 0.4, Δz = 0.5, ψ_sat = -0.1, ρ_liq = 1000. Surfaces the supercooled-water response surface — w_liq_max increasing toward the upper-right (warmer T_i and more clay-like B_i) and shrinking toward the lower-left (colder T_i and more sand-like B_i).",
      "initial_state": {},
      "parameters": {
        "θ_sat": 0.4,
        "Δz": 0.5,
        "ψ_sat": -0.1,
        "ρ_liq": 1000.0
      },
      "time_span": {
        "start": 0.0,
        "end": 1.0
      },
      "parameter_sweep": {
        "type": "cartesian",
        "dimensions": [
          {
            "parameter": "T_i",
            "range": {
              "start": 250.0,
              "stop": 273.0,
              "count": 12,
              "scale": "linear"
            }
          },
          {
            "parameter": "B_i",
            "range": {
              "start": 3.0,
              "stop": 12.0,
              "count": 10,
              "scale": "linear"
            }
          }
        ]
      },
      "plots": [
        {
          "id": "w_liq_max_surface",
          "type": "heatmap",
          "description": "w_liq_max as a function of T_i (cold branch) and B_i: the supercooled-water response surface per Eq. 4.58, monotone in both axes.",
          "x": {
            "variable": "T_i",
            "label": "T_i (K)"
          },
          "y": {
            "variable": "B_i",
            "label": "B_i (dimensionless)"
          },
          "value": {
            "variable": "w_liq_max",
            "reduce": "final"
          }
        }
      ]
    }
  ]
}