Extend climate and land-use drivers back to 1750 for WHEP's legacy processes (SOC, N, future P, LPJmL)

## Summary

WHEP's gridded climate drivers start in **1901** (CRU) and its land-use build
starts in **1850**. The canonical WHEP period is **1850–2023**. Every process in
WHEP with long memory — soil organic carbon, soil organic nitrogen, legacy
phosphorus, hydrology, and everything LPJmL simulates — therefore begins its
history at, or after, the first year we want to report on. The state entering
1850 is not simulated; it is assumed.

This issue asks for **climate and land-use drivers extended back to 1750**,
giving roughly a century of real, transient history before the canonical period
opens. The reconstruction should use the **best available source for each period
and variable** (a reliable series starting in 1820 should be used from 1820,
with a weaker source only below that), and must preserve **time-series
coherence** by splicing **anomalies relative to a common normal period**, never
raw values from mixed sources.

This is driver-data work, not model work. LPJmL is one consumer among several.

## Why: legacy processes cannot start at the year we want to report

**Soil organic carbon.** [`build_carbon_balance()`](R/carbon_balance.R)
initialises every cell by running the turnover model to **equilibrium under the
earliest year's per-land-use carbon inputs**, then weighting those equilibria by
the earliest year's land-use fractions
([R/carbon_balance.R:428-441](R/carbon_balance.R#L428-L441)). Starting the march
in 1850 therefore asserts that global soils were *at steady state in 1850* —
which is false for exactly the regions that matter most. European, Chinese and
Indian croplands had been under cultivation for centuries by 1850 and were
mid-decline, not at equilibrium. The resulting initialisation error propagates
through the entire canonical period as a spurious ΔSOC trend.

**Soil organic nitrogen.** The N balance derives ΔSON from ΔSOC, so it inherits
the same initialisation artefact directly.

**Legacy phosphorus (future work).** Soil P accumulates over more than a century
and its legacy pool is the whole point of modelling it. A P module that starts in
1850 with an assumed equilibrium pool cannot answer the question it exists to
answer. For P a pre-1850 run-in is a hard requirement, not a refinement — which
is why the driver data should be built once, now, rather than per module.

**LPJmL.** Every WHEP-driven run starts in **1901**, because that is where CRU
begins ([run_lpjml.R:26](inst/scripts/run_lpjml.R#L26)). NPP, hydrology/drainage
and grass productivity for 1850–1900 are thus spin-up output rather than
simulation, and 1901–1950 still carries initialisation drift.

A shared 1750-onward driver set fixes all four at once, and the 1750–1850 window
contains real forced signals that a recycled climatology erases entirely: the
Dalton Minimum, the 1809 eruption, Tambora (1815) and the 1816 "year without a
summer".

## Current state: where each driver stops today

| Driver | Source today | Starts | Consumers |
|---|---|---|---|
| Temperature, precipitation (monthly) | CRU TS ([R/cru_climate.R](R/cru_climate.R)) | 1901 | SOC turnover modifiers, LPJmL, future P |
| Cloud cover, wet days | CRU TS | 1901 | LPJmL |
| Wind | GSWP3-W5E5 | 1901 | LPJmL |
| Land-use classes | WHEP build (LUH2 v2h 850–2015 available, [R/luh2_landuse.R](R/luh2_landuse.R)) | 1850 | SOC march + LUC transfer, LPJmL |
| Carbon inputs / management | WHEP build | 1850 | SOC |
| CO2 | `historical_CO2_annual_1765_2018.txt` | 1765 | LPJmL |
| N deposition (NHx/NOy) | HaNi | 1860 | N balance, LPJmL |
| Soil, lakes/rivers | static | n/a | — |

The concrete LPJmL file names and their year spans are in
[run_lpjml.R:66-75](inst/scripts/run_lpjml.R#L66-L75). Note that
[R/cru_climate.R](R/cru_climate.R) reads CRU TS 4.09 (1901–2024) while the LPJmL
run is still pinned to CRU TS 3.10 (1901–2009); both share the same 1901 wall,
so this issue is orthogonal to which CRU release is used.

## Design principle: anomalies over normals, never raw splicing

For every variable, the reconstruction must be built as

> `value(t) = CRU_normal(cell, month) (+ or *) anomaly_source(t, cell, month)`

where `anomaly_source` is the older dataset expressed as a departure from **its
own** normal over a common overlap window, and `CRU_normal` is the CRU
climatology over that same window. Additive for temperature, multiplicative
(with physical bounds) for precipitation, cloud, wet days and wind.

This buys three things that matter more than the absolute skill of any single
source:

1. **No step at the splice.** By construction the reconstruction converges to
CRU at the overlap, so 1900→1901 shows no discontinuity in the mean.
2. **CRU's spatial detail is preserved** even where the older source is coarse
(ModE-RA and EKF400 are ~2 degrees, 20CRv3 ~75 km, all much coarser than the
0.5 degree WHEP grid).
3. **Sources can change between periods without a jump.** If one source covers
1806–1900 and a weaker one 1750–1805, both are expressed as anomalies on the
same CRU normal, so the transition between them is smooth as long as their
own overlap is used to cross-calibrate.

The overlap/normal window is itself an open decision (see below).

## Land use is the delicate case: no 1850 jump

WHEP's own land-use build starts in 1850. LUH2 v2h covers 850–2015 and is
already read by [R/luh2_landuse.R](R/luh2_landuse.R), so the raw data for
1750–1850 exists. **Do not splice raw LUH2 onto WHEP at 1850** — the two differ
systematically in cropland and grassland extent, and a raw splice would put an
artificial step in the land-use series at the first year of the canonical
period, which `build_carbon_balance()` would read as a real land-use change flux
and convert into a spurious LUC carbon transfer.

The same anomaly-over-normal principle applies: take LUH2's **relative change**
over 1750–1850 and apply it backwards from the WHEP 1850 state, per cell and per
land-use class, so the reconstructed 1850 land use equals WHEP's 1850 land use
by construction. Class-area consistency within a cell (cropland + grassland +
natural + urban) must be preserved after the rescaling.

Carbon inputs and management for 1750–1850 need the same treatment and are
likely the weakest link: pre-1850 yields, residue returns and manure are far
less constrained than land-use extent. Say explicitly what is assumed.

## Open methodological decisions (needs-expert)

These embed scientific choices and must not be settled on green CI alone.

1. **What the equilibrium initialisation should be conditioned on once the
transient starts in 1750.** `build_carbon_balance()` already accepts an
`equilibrium_climate` pre-industrial normal separate from the forward march
([R/carbon_balance.R:444-449](R/carbon_balance.R#L444-L449)), and LPJmL's
5000-year spin-up currently recycles detrended 1901–1930 CRU. Two options:
- Detrended **1750–1780 reconstruction** — pre-industrial-consistent: pools
equilibrate to the same climate the transient then begins from. Costs: it
leans on the reconstruction at its oldest and weakest end.
- Keep **1901–1930** — the status quo and the LPJmL/TRENDY convention. Costs:
pools equilibrate to 20th-century climate, then the transient starts 150
years earlier.

Decide at implementation time, ideally with a sensitivity run.
2. **Normal / overlap window** for the anomaly splice (e.g. 1901–1930 vs
1961–1990). The earlier window is closer to the reconstruction period and
overlaps more sources; the later is the WMO convention and better observed.
3. **Variables with no pre-1901 observational basis** (cloud cover, wet days,
wind). Options: take them from a reanalysis that carries them (20CRv3),
derive them from reconstructed precipitation and temperature via the
relationships CRU itself exhibits, or hold them at a detrended recycled
climatology while temperature and precipitation vary. State the choice per
variable; do not silently fall back.
4. **Pre-1860 N deposition.** HaNi starts 1860. Scaling its 1860 field toward a
natural pre-industrial background over 1750–1860 is one option; a
chemistry-based historical deposition product reaching further back is
another. Either way the 1860 value must match HaNi exactly at the join.

## Recommended sources to evaluate (not prescriptive)

The implementer selects the final sources and methods. These are verified
candidates worth checking first, per the reuse-first rule.

**Climate (the binding constraint)**

- **ModE-RA** — global monthly paleo-reanalysis, **1421–2008**, offline data
assimilation blending an AGCM ensemble with proxy, documentary and early
instrumental data. Covers the whole target window in one product, which is the
strongest argument for it. Valler et al. (2024), *Scientific Data*,
[doi:10.1038/s41597-023-02733-8](https://doi.org/10.1038/s41597-023-02733-8).
- **20CRv3** — NOAA-CIRES-DOE 20th Century Reanalysis v3, **1806–2015**
(1806–1835 flagged experimental), 3-hourly, T254 Gaussian (~75 km), 80
members, assimilating surface pressure only. Highest resolution and the only
candidate natively carrying cloud, wind and humidity, but it does not reach
1750.
- **EKF400v2** — monthly paleo-reanalysis **1603–2003** at 2 degrees, Kalman
filtering of early instrumental, documentary and tree-ring data into an AGCM
ensemble. Franke et al. (2020), WDCC,
[doi:10.26050/WDCC/EKF400_v2.0](https://doi.org/10.26050/WDCC/EKF400_v2.0).
- **Berkeley Earth land TAVG** — monthly gridded 1-degree temperature
**anomalies from 1750**, referenced to 1951–1980. Already in anomaly form,
which fits the splice design directly, and the longest instrumental
temperature field available. Caveat: CC BY-NC 4.0 (non-commercial) — check
licence compatibility before depending on it.
- **GPCC Full Data Monthly v2022** — gauge-based land precipitation,
**1891–2020**, available at 0.5 degrees. Extends only 10 years beyond CRU, but
those are real gauge observations and make an excellent independent check on
reconstructed precipitation skill at the splice. Schneider et al. (2022),
[doi:10.5676/DWD_GPCC/FD_M_V2022_050](https://doi.org/10.5676/DWD_GPCC/FD_M_V2022_050).

**Land use and other drivers**

- **Land use (pre-1850)** — LUH2 v2h (Hurtt et al. 2020, *GMD*), 850–2015 at
0.25 degrees, already wired in [R/luh2_landuse.R](R/luh2_landuse.R). Used as a
*relative change* driver anchored on WHEP's 1850 state, per the section above.
- **CO2 (pre-1765)** — Meinshausen et al. (2017), *GMD* 10, 2057–2116,
[doi:10.5194/gmd-10-2057-2017](https://doi.org/10.5194/gmd-10-2057-2017), gives
historical GHG concentrations from year 0 to 2014, so 1750–1764 is a direct
read. Verify agreement with the current `historical_CO2_annual_1765_2018.txt`
over 1765–2014 before splicing.
- **N deposition (pre-1860)** — HaNi (Tian et al. 2022, *ESSD* 14, 4551–4568,
[doi:10.5194/essd-14-4551-2022](https://doi.org/10.5194/essd-14-4551-2022);
data [doi:10.1594/PANGAEA.942069](https://doi.org/10.1594/PANGAEA.942069)) is
the current source and starts 1860; see decision 4 above.

## Scope

A full builder in WHEP: readers for the selected pre-1901 sources, an
anomaly-splicing builder, and writers emitting the extended drivers in **both**
forms WHEP needs — the tidy tibble contract the SOC/N modules consume
(`get_soc_climate_drivers()`, the `land_use` input of `build_carbon_balance()`)
and LPJmL-ready NetCDF spliced onto the existing CRU/GSWP3 files at 1901. One
driver set, several consumers; a future P module should be able to use it
unchanged.

Following the multi-method convention in `CLAUDE.md`, the builder should expose
a `method =` / `source =` argument selecting among candidate reconstructions
rather than hardcoding one, default to the most defensible, and **record the
contributing source per variable-period in an output column and in the NetCDF
attributes**. Sources are alternatives, never silent fallbacks.

## Acceptance criteria

- [ ] Climate and land-use drivers are available for **1750 onward** through the
WHEP reader API, and `build_carbon_balance()` can march from 1750.
- [ ] `run_lpjml()` can be configured with `simulation_start_year = 1750` and
completes, with `.check_climate_coverage()` passing for the full period.
- [ ] **No splice discontinuity**: the 30-year running mean of each variable
shows no step at 1900→1901 beyond the interannual noise level, verified
globally and for a set of contrasting regions.
- [ ] **Physical bounds hold** for every reconstructed cell-month: precipitation
>= 0, cloud within 0–100%, wet days <= days in month, wind > 0, temperature
within a plausible range.
- [ ] **Known events are reproduced**: the expected cooling signal around 1816
and the Dalton Minimum are visible, i.e. the run-in carries real forced
variability and is not a recycled climatology.
- [ ] **Land use is continuous at 1850**: the reconstructed 1850 state equals
WHEP's 1850 state per cell and class, and `build_carbon_balance()` reports
no anomalous LUC transfer flux in 1850.
- [ ] CO2 at 1765 and N deposition at 1860 match their existing sources exactly
at the join.
- [ ] A 1750-start vs 1850-start SOC run is compared, quantifying how much of
the 1850–1900 ΔSOC signal was initialisation artefact.
- [ ] Source directories resolve from environment variables (never hardcoded
absolute paths), consistent with the existing raster readers.
- [ ] Tests with small fixtures covering the splice arithmetic, the physical
bounds and the land-use anchoring; `@examples` use the `example = TRUE`
fixture pattern.
- [ ] Chosen sources, normal window and equilibrium decision documented with
verified citations.

## Related

- #340 — the other end of the same time axis: LPJmL output effectively ends
2009, so covering the full canonical 1850–2023 period needs both this issue
and a usable post-2009 run.
- #351, #358 — the SOC modules' climate and clay readers are unwired stubs;
whatever reader API this issue lands should serve them.
- #352 — closed-form equilibrium spin-up; interacts with decision 1 above.

cc @lbm364dl @t-sum97 @Jauviltra