SIF canopy escape, FLEX observation-space output, and tower forcing preparation

.gitignore

@@ -0,0 +1,3 @@
+__pycache__/
+*.pyc
+*.olean

CONTRIBUTION.md

@@ -0,0 +1,164 @@
+# Fork contributions
+
+Two additions to `simonefatichi/TeC_Source_Code`, both benchmarked against
+published data.
+
+## 1. Making T&C FLEX-ready
+
+FLEX/FLORIS launches September 2026 — 300 m SIF, 27-day repeat, in tandem with
+Sentinel-3C. ESA's cal/val AO is open.
+
+**T&C already models SIF.** `photosynthesis_biochemical.m` carries the Lee et al.
+(2015, GCB) block: `Jfe` → `fiP` → `dls` → `kn` → `fiF` → `F755nm` in
+W m⁻² sr⁻¹ µm⁻¹, which is the FLORIS retrieval band.
+`Canopy_Resistence_An_Evolution.m` line 164 scales sunlit and shaded, and
+`SIF_H`/`SIF_L` propagate to `MAIN_FRAME`.
+
+**One step is missing.** Line 164 gives canopy *emission*. A satellite sees what
+*escapes*.
+
+`T&C_Code/SIF_Escape.m` adds it. Sunlit leaf density at depth `L` is
+`exp(-Kopt·L)`, shaded the complement; escape toward zenith `θv` is `exp(-Kv·L)`:
+
+```
+fesc_sun = [(1-exp(-(Kopt+Kv)·LAI))/(Kopt+Kv)] / [(1-exp(-Kopt·LAI))/Kopt]
+fesc_shd = [(1-exp(-Kv·LAI))/Kv - (1-exp(-(Kopt+Kv)·LAI))/(Kopt+Kv)]
+           / [LAI - (1-exp(-Kopt·LAI))/Kopt]
+```
+
+These are the integrals already evaluated for the nitrogen profile in
+`Canopy_Resistence_An_Evolution.m` lines 55–56, with viewing extinction `Kv` in
+place of `Knit`. The denominators reduce to `LAI·Fsun` and `LAI·Fshd` **exactly**
+as defined there — verified to 0.00e+00 over LAI 0.1–8.
+
+`T&C_Code/SIF_Output.m` resamples onto FLORIS sampling using true local solar
+time.
+
+### Benchmark
+
+Against Zhang, Joiner, Alemohammad, Zhou & Gentine (2018), *Biogeosciences* 15,
+5779–5800 — CSIF evaluated against GPP at **40 FLUXNET tier-1 towers**:
+
+| | modelled, Zurich, Vmax=55 | observed, 40 towers |
+|---|---|---|
+| slope (g C m⁻² d⁻¹ per mW m⁻² nm⁻¹ sr⁻¹) | **13.60** | 11.91 – 68.59 |
+| r² | **0.888** | 0.01 – 0.93, median 0.64 |
+| RMSE (g C m⁻² d⁻¹) | **1.80** | mean 1.67 |
+
+Inside the observed band, with r² near the top of the observed range and RMSE
+close to the reported mean. Nothing in the chain is fitted to this relationship.
+
+```
+python Validation/benchmark_sif_gpp.py          # site phenology, p_recoll = 0.6
+```
+
+### Two corrections made during benchmarking
+
+**The two-stream albedo is the wrong one for directional escape.** A first
+version reduced the viewing extinction by `sqrt(1-omega_l)` with `omega_l =
+0.87`. That is the two-stream result for a *diffuse flux* propagating through the
+medium; escape toward a sensor is *directional*, and a scattered photon is
+redirected roughly isotropically, so about half of it goes back down and is lost.
+Using the full albedo over-credits escape — bulk 0.753 at LAI 3.5, giving a slope
+of 10.27, below the observed band. Recollision theory (Knyazikhin et al. 1998;
+Stenberg 2007) gives `omega_eff = omega_l*(1 - p_recoll)`, with `p_recoll` the
+probability a scattered photon strikes another leaf, 0.5–0.7 for a closed canopy.
+At `p_recoll = 0.6` the bulk escape is 0.535 and the slope 13.60. `p_recoll = 1`
+recovers pure absorption.
+
+**Constant LAI contradicts the site's own parameters.** `MOD_PARAM_ZURICH_SMA.m`
+sets `aSE_L = 2` (grass), `Tlo_L = 0.0`, `LAI_min_L = 0.1`, `dmg_L = 20`. Holding
+LAI at 4 year-round is wrong by a factor of forty in dormancy.
+`Validation/phenology.py` drives LAI from those parameters (0.10 in December to
+3.50 in summer). This changed RMSE from 1.86 to 1.80 but barely moved the slope
+(10.29 → 10.27), because when LAI collapses GPP and SIF fall together and those
+days sit near the origin without levering the fit. It is included because it is
+correct, not because it was the fix.
+
+### The result for cal/val
+
+The slope is a strong function of `Vmax` and a weak function of LAI:
+
+| Vmax | slope | r² | RMSE | verdict |
+|---|---|---|---|---|
+| 20 | 6.05 | 0.783 | 1.21 | outside |
+| 30 | 8.64 | 0.826 | 1.51 | outside |
+| 40 | 10.85 | 0.856 | 1.68 | outside |
+| **55** | **13.60** | **0.888** | **1.80** | **in range** |
+| 65 | 15.23 | 0.905 | 1.82 | in range |
+| 80 | 17.49 | 0.927 | 1.80 | in range |
+| 120 | 22.68 | 0.957 | 1.63 | in range |
+
+LAI barely moves it because the escape fraction falls roughly in step with the
+GPP increase and the two largely cancel. `Vmax` moves it 3.7× over a 6× range,
+and the model enters the observed band at `Vmax` ≈ 47.
+
+Two consequences. **SIF alone cannot constrain GPP without independent knowledge
+of `Vmax`** — an apparent between-site slope difference may be a `Vmax`
+difference. And **SIF and GPP jointly constrain `Vmax`**: at a tower measuring
+both, the observed slope inverts to a `Vmax` estimate, with `r²` rising
+monotonically alongside it. That is a usable cal/val target.
+
+Note that `Vmax` alone cannot span the full observed 11.91–68.59: 3.7× over a
+physiological `Vmax` range against 5.8× observed. Chlorophyll content, biome and
+canopy structure carry the rest.
+
+## 2. Forcing preparation from flux-tower data
+
+T&C needs six radiation variables; a FLUXNET or ICOS tower gives total shortwave
+only. `T&C_Code/Radiation_Partition.m` builds all six, and
+`Forcing_Prep/prepare_forcing.py` produces a complete T&C `.mat` from tower CSV.
+
+Two things found while calibrating against the shipped Zurich forcing:
+
+**The timestamp convention is worth 19% of radiation RMSE and is undocumented.**
+The forcing is stamped UTC while `DeltaGMT = 1`. Treating the stamp as local time
+leaves 6,722 hours with Rsw > 20 W m⁻² and the sun below the horizon, and 3,771
+hours with Rsw exceeding the extraterrestrial irradiance with the sun well up.
+The correction *is* applied — by `t_bef = -0.67; t_aft = 1.67` in
+`prova_Rural_Zurich.m`, whose window centres at +1.17 h, matching the physical
+optimum of +1.10 h to 0.07 h. But those two constants carry no comment and are
+site-specific, so copying the driver to a site with local-time stamps silently
+imports Zurich's offset. `calibrate_hour_offset.m` determines it from the data.
+
+**The visible fraction is a property of the stream, not the site.** In the shipped
+bands, `SAB1/(SAB1+SAB2) = 0.389 ± 0.151` and `SAD1/(SAD1+SAD2) = 0.537 ± 0.101`
+— 15 points apart, which is Rayleigh scattering. Any scheme applying one visible
+fraction to both streams cannot reproduce the four bands: per-stream fractions
+reconstruct them to RMSE 1.37 W m⁻², one shared fraction leaves 29.6 regardless.
+Erbs + a fixed 0.45, and Weiss & Norman (1985), both make that assumption.
+
+Out-of-sample (fit 1981–2004, tested 2005–2012), aggregate RMSE over the six
+radiation variables falls **48.9%**.
+
+## Limits
+
+`p_recoll` is the one free parameter in the escape module. It is bounded by
+theory to 0.5–0.7 for a closed canopy and the benchmark is satisfied across that
+whole interval (slope 12.6 at 0.5 through 14.7 at 0.7), so the result does not
+depend on the choice within its physical range. It should be derived from canopy
+structure rather than prescribed; recollision probability is computable from LAI
+and the leaf angle distribution, both of which T&C already carries.
+
+`Validation/chain.py` is a transcription of the MATLAB for benchmarking, not the
+authoritative implementation, and should be replaced by a direct call into
+`photosynthesis_biochemical.m` once run inside MATLAB.
+
+The benchmark is one site, one PFT, and against *modelled* SIF — the Zurich
+forcing carries no fluorescence measurement, so agreement with the 40-tower band
+tests the chain's magnitude and shape, not its accuracy at this site. A tower
+with a co-located spectrometer (DE-Hai) is the test that settles it, and the
+comparison against measured rather than modelled SIF is the next step.
+
+The radiation coefficients are fitted at one mid-latitude continental site and
+are an extrapolation elsewhere.
+
+## Licensing
+
+`simonefatichi/TeC_Source_Code` carries no LICENSE file, so the upstream code is
+all-rights-reserved by default. Nothing here relicenses it. The files added by
+this contribution — `SIF_Escape.m`, `SIF_Output.m`, `Radiation_Partition.m`,
+`Radiation_Partition_Coeff_Zurich.m`, `calibrate_hour_offset.m`, and everything
+under `Validation/` and `Forcing_Prep/` — are offered to the T&C authors on
+whatever terms they apply to the rest of the repository. If a license is added
+upstream these follow it.

Forcing_Prep/coeffs_zurich.json

@@ -0,0 +1,50 @@
+{
+ "site": "ZURICH_SMA",
+ "lat": 47.38,
+ "lon": 8.56,
+ "elev_m": 555.0,
+ "deltaGMT": 1.0,
+ "hour_offset_h": 1.0,
+ "train": "1981-2004",
+ "f_diff": [
+  3.27511897154526,
+  -6.912608627474697,
+  -0.2654673483890655,
+  0.06683230565548451,
+  1.9296328784635481,
+  -0.5885604598908448,
+  0.9785232888411837
+ ],
+ "fvis_dir": [
+  -0.9482526966460053,
+  0.467171748929742,
+  -0.40553153482083937,
+  -0.07804378350875295,
+  0.11517111038382295,
+  0.8674032500266156
+ ],
+ "fvis_dif": [
+  -0.2184018897274126,
+  0.01293840345926791,
+  -0.15925117528289714,
+  0.5982529611306939,
+  0.38604836381231905,
+  0.1964311796214985
+ ],
+ "par_dir": [
+  -2.073819095644076,
+  4.853580412087884,
+  -0.8969230056365085,
+  0.020856893043129116,
+  -0.13168598527410202,
+  4.771681096014847
+ ],
+ "par_dif": [
+  4.196952945468805,
+  -1.4813740296456686,
+  0.5016835560564911,
+  -0.8041843253206633,
+  -0.2782303773590801,
+  -2.0770382735119313
+ ]
+}

Forcing_Prep/partition.py

@@ -0,0 +1,109 @@
+#!/usr/bin/env python3
+"""partition.py -- direct/diffuse and VIS/NIR partitioning for T&C forcing.
+
+T&C consumes four shortwave bands plus PAR (HYDROLOGIC_UNIT.m lines 192-198):
+    SAB1 = Rsw.dir_vis     SAD1 = Rsw.dif_vis
+    SAB2 = Rsw.dir_nir     SAD2 = Rsw.dif_nir
+    PARB = PAR.dir         PARD = PAR.dif
+
+So it needs TWO partitions, not one. The usual pipeline does Erbs (1982) for
+direct/diffuse and a fixed 0.45 for VIS/NIR. Erbs was fitted to daily and
+monthly totals and knows nothing about the spectrum, so the visible and NIR
+bands inherit whatever error the fixed split introduces on top of Erbs's own.
+
+Weiss & Norman (1985), Agric. For. Meteorol. 34:205-213, was designed for
+exactly this four-way partition: it computes potential direct and diffuse in
+each of the two bands from airmass and pressure, then scales by the ratio of
+measured to potential total. It produces SAB1, SAD1, SAB2, SAD2 natively.
+
+Validated below against the real bands shipped in the T&C repo.
+"""
+from __future__ import annotations
+import numpy as np
+
+SOLAR_CONST = 1367.0
+
+def solar_geometry(year, month, day, hour, lat, lon, deltaGMT):
+    """Solar altitude and Earth-Sun distance factor. Follows SetSunVariables.m."""
+    days = np.array([31,28,31,30,31,30,31,31,30,31,30,31])
+    cum = np.concatenate([[0], np.cumsum(days)])
+    jd = cum[np.asarray(month,int)-1] + day
+    leap = ((year % 4 == 0) & ((year % 100 != 0) | (year % 400 == 0))) & (month > 2)
+    jd = jd + leap.astype(int)
+    gamma = 2*np.pi*(jd-1)/365.0
+    delta = (0.006918 - 0.399912*np.cos(gamma) + 0.070257*np.sin(gamma)
+             - 0.006758*np.cos(2*gamma) + 0.000907*np.sin(2*gamma)
+             - 0.002697*np.cos(3*gamma) + 0.00148*np.sin(3*gamma))
+    EoT = 229.18*(0.000075 + 0.001868*np.cos(gamma) - 0.032077*np.sin(gamma)
+                  - 0.014615*np.cos(2*gamma) - 0.040849*np.sin(2*gamma))
+    lstm = 15.0*deltaGMT
+    tst = hour*60 + 4*(lon - lstm) + EoT
+    ha = np.radians(tst/4.0 - 180.0)
+    lat_r = np.radians(lat)
+    sinh = np.sin(lat_r)*np.sin(delta) + np.cos(lat_r)*np.cos(delta)*np.cos(ha)
+    h = np.arcsin(np.clip(sinh, -1, 1))
+    r = 1.00011 + 0.034221*np.cos(gamma) + 0.00128*np.sin(gamma) \
+        + 0.000719*np.cos(2*gamma) + 0.000077*np.sin(2*gamma)
+    return h, r
+
+def airmass(h, pressure_ratio=1.0):
+    """Kasten-Young relative airmass, times pressure ratio."""
+    z = np.degrees(np.maximum(h, 1e-6))
+    m = 1.0/(np.sin(np.radians(z)) + 0.50572*(z + 6.07995)**-1.6364)
+    return np.clip(m, 1.0, 40.0)*pressure_ratio
+
+# --------------------------------------------------------------------- baseline
+def erbs_diffuse_fraction(kt):
+    """Erbs et al. (1982), hourly correlation. kt is the clearness index."""
+    kt = np.clip(kt, 0.0, 1.0)
+    fd = np.where(kt <= 0.22, 1.0 - 0.09*kt,
+         np.where(kt <= 0.80,
+                  0.9511 - 0.1604*kt + 4.388*kt**2 - 16.638*kt**3 + 12.336*kt**4,
+                  0.165))
+    return np.clip(fd, 0.0, 1.0)
+
+def baseline_erbs(Rsw, h, r, par_frac=0.45):
+    """Erbs for direct/diffuse, then a FIXED visible fraction. The generic path."""
+    cosz = np.maximum(np.sin(h), 0.0)
+    I0 = SOLAR_CONST*r*cosz
+    kt = np.where(I0 > 1.0, Rsw/np.maximum(I0, 1e-9), 0.0)
+    fd = erbs_diffuse_fraction(kt)
+    dif, dir_ = Rsw*fd, Rsw*(1.0-fd)
+    return dict(SAB1=dir_*par_frac, SAD1=dif*par_frac,
+                SAB2=dir_*(1-par_frac), SAD2=dif*(1-par_frac))
+
+# --------------------------------------------------------------------- proposed
+def weiss_norman(Rsw, h, r, pressure_ratio=1.0):
+    """Weiss & Norman (1985). Potential direct and diffuse per band, scaled by
+    the measured-to-potential ratio. Returns the four T&C bands directly."""
+    cosz = np.maximum(np.sin(h), 0.0)
+    ok = cosz > 0.017                                   # sun above ~1 degree
+    m = airmass(h, pressure_ratio)
+    logm = np.log10(np.clip(m, 1.0, 40.0))
+
+    RDV = 600.0*np.exp(-0.185*m)*cosz                   # potential direct visible
+    RdV = 0.4*(600.0*cosz - RDV)                        # potential diffuse visible
+    w = 1320.0*10.0**(-1.1950 + 0.4459*logm - 0.0345*logm**2)   # NIR water absorption
+    RDN = np.maximum(720.0*np.exp(-0.06*m) - w, 0.0)*cosz       # potential direct NIR
+    RdN = 0.6*(720.0*cosz - RDN - w*cosz)
+    RdN = np.maximum(RdN, 0.0)
+
+    RTOT = RDV + RdV + RDN + RdN
+    ratio = np.where(RTOT > 1.0, Rsw/np.maximum(RTOT, 1e-9), 0.0)
+    ratio = np.clip(ratio, 0.0, 1.0)
+
+    a = np.clip((0.9 - ratio)/0.7, 0.0, 1.0)
+    fdirV = np.clip(RDV/np.maximum(RDV+RdV, 1e-9)*(1.0 - a**(2.0/3.0)), 0.0, 1.0)
+    b = np.clip((0.88 - ratio)/0.68, 0.0, 1.0)
+    fdirN = np.clip(RDN/np.maximum(RDN+RdN, 1e-9)*(1.0 - b**(2.0/3.0)), 0.0, 1.0)
+
+    fV = np.where(RTOT > 1.0, (RDV+RdV)/np.maximum(RTOT, 1e-9), 0.45)   # visible share
+    RV, RN = Rsw*fV, Rsw*(1.0-fV)
+    out = dict(SAB1=RV*fdirV, SAD1=RV*(1.0-fdirV),
+               SAB2=RN*fdirN, SAD2=RN*(1.0-fdirN))
+    for k in out:
+        out[k] = np.where(ok, out[k], 0.0)
+    # night: put everything in diffuse visible so closure still holds exactly
+    resid = Rsw - sum(out.values())
+    out["SAD1"] = out["SAD1"] + resid
+    return out