diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3edefc9..b24a12c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,12 @@
 # Changelog
 
+## Repeated calendar-block spin-up — 2026-10-20
+
+- Add complete, consecutive `spinup.source_years` blocks to both solver backends, calibration and native workflows, retaining leap days and the original output calendar.
+- Count repetitions including the retained source block and exclude that block from scoring. Keep memory independent of repetition count by discarding each column's extra warm-up trajectories.
+- Retain mandatory convergence gates by default; support explicitly labelled fixed-duration experiments that preserve failed diagnostics without claiming convergence.
+- Validate compact histories against expanded histories and verify TOML/API integration; 480 offline tests passed in Gadi job 181730321.
+
 ## Daily two-parameter SWB — 2026-10-18
 
 - Use one conservative update per daily forcing record; retain donor/receiver limits, physical bounds and separate actual/unmet E/T budgets.
diff --git a/README.md b/README.md
index 08920dc..57fc703 100755
--- a/README.md
+++ b/README.md
@@ -61,7 +61,9 @@ The five stages share one implementation in `pysweb.workflow`:
 | `calibrate` | Fit two domain parameters to reference surface soil moisture | Enabled explicitly; local reference or configured Earth Engine asset |
 | `run` | Simulate independent soil columns | Prepared forcing and hydraulic properties |
 
-Select stages with `--stages preprocess,calibrate,run`. The default workflow preserves source download grids, fits on native GSSM cells, and writes both final products on one selected Landsat 30 m UTM grid. It runs only Landsat 8/9 by default, for dates from 2013 onwards. Earlier missions require explicit `run.sensors = ["LT05", "LE07", "LC08", "LC09"]`; no second sensor pathway is launched automatically. For an 8/9 run beginning in 2013, calibration starts on 2014-01-01 by default, reserving the available 2013 record for initialisation. Set `calibration.start` to override this choice. For both calibration and normal runs, the SWB makes one conservative update per day, starts at physical field capacity and carries the available 2013 warm-up state into 2014. Resolved plans and final products identify the warm-up; NetCDF `is_warmup` flags those dates. No years are repeated unless an optional `[spinup]` table is explicitly supplied. Verify forcing coverage; this initialization assumes field capacity rather than claiming complete independence from the starting profile. Calibration is skipped when `calibration.enabled = false`. Its scoring window lies within the continuous simulation window. Set `run.holdout_start` to protect an independent test period. See [native grids and output contracts](docs/native-grids.md). Optional [configured acceptance protocols](docs/configured-protocols.md) integrate validated L7 donors, development-only spin-up, paired support, benchmark timing and site-specific evaluation. See [the workflow guide](docs/workflows.md) for offline operation, HPC execution and individual APIs.
+Select stages with `--stages preprocess,calibrate,run`. The default workflow preserves source download grids, fits on native GSSM cells, and writes both final products on one selected Landsat 30 m UTM grid. It runs only Landsat 8/9 by default, for dates from 2013 onwards. Earlier missions require explicit `run.sensors = ["LT05", "LE07", "LC08", "LC09"]`; no second sensor pathway is launched automatically. For an 8/9 run beginning in 2013, calibration starts on 2014-01-01 by default, reserving the available 2013 record for initialisation. Set `calibration.start` to override this choice. For both calibration and normal runs, the SWB makes one conservative update per day, starts at physical field capacity and carries the available 2013 warm-up state into 2014. Resolved plans and final products identify the warm-up; NetCDF `is_warmup` flags those dates. Repeated spin-up is disabled by default; set `[spinup] enabled = true` to enable it, or `enabled = false` to retain the settings without repetition. Verify forcing coverage; this initialization assumes field capacity rather than claiming complete independence from the starting profile. Calibration is skipped when `calibration.enabled = false`. Its scoring window lies within the continuous simulation window. Set `run.holdout_start` to protect an independent test period. See [native grids and output contracts](docs/native-grids.md). Optional [configured acceptance protocols](docs/configured-protocols.md) integrate validated L7 donors, development-only spin-up, paired support, benchmark timing and site-specific evaluation. See [the workflow guide](docs/workflows.md) for offline operation, HPC execution and individual APIs.
+
+For repeated multi-year initialization, see [calendar-block spin-up](docs/multiyear-spinup.md). A complete 2014–2015 block can be repeated five or ten times without expanding forcing grids or spatial output calendars. Repetition counts include the retained source block; scoring starts afterward. Convergence is checked separately from numerical conservation.
 
 Every configured run writes `resolved-config.json`. Model output manifests record input SHA-256 checksums, configuration, code identity and output checksums. `skip_existing` validates the model output against this identity before reuse. Unmanifested older outputs require regeneration. Generated data and environments stay outside version control.
 
@@ -124,4 +126,4 @@ The code is distributed under the [MIT licence](LICENSE). Dataset licences remai
 
 Component outputs use `E_s` (soil evaporation), `E_t` (transpiration) and `f_t` (transpiration fraction). See the [notation migration and legacy-input policy](docs/migration.md#component-notation-migration-11-october-2026).
 
-The [conservative two-parameter SWB formulation](docs/swb-formulation.md) calibrates diffusivity and root beta, uses fixed versioned GSSM references and records separate actual/unmet E/T with conservative water budgets. Existing four-parameter calibrations require recalibration.
+The [conservative two-parameter SWB formulation](docs/swb-formulation.md) calibrates diffusivity and root beta, uses fixed versioned GSSM references and records separate actual/unmet E/T with conservative water budgets. V3 applies the GSSM/WP dry ratio to extraction thresholds in 0–15 cm, a half adjustment in 15–30 cm, and hydraulic WP below 30 cm. Continuous three-day relaxation and physical zero-to-porosity bounds are retained. Existing v2 and four-parameter calibrations require recalibration; immutable GSSM reference artifacts remain compatible.
diff --git a/docs/configured-protocols.md b/docs/configured-protocols.md
index 0ea532a..3049ed6 100644
--- a/docs/configured-protocols.md
+++ b/docs/configured-protocols.md
@@ -22,12 +22,17 @@ The padding table contains `primary_bounds`, `comparison_bounds`, their `crs`, `
 
 ## Paired support and initialization
 
-The usual Landsat 8/9 workflow runs alone. `forcing.peer_runs` is optional and should be set only for an explicitly requested comparison. For runs beginning in 2013, the default calibration scoring start is 2014-01-01, so available 2013 forcing supplies initialisation. Both calibration and normal runs (including `calibration.enabled = false`) start at physical field capacity, run the available 2013 record once and carry its state into 2014 without resetting. Resolved plans record the initialization and analysis windows; final SWB products preserve the actual calendar and mark warm-up days with `is_warmup`. Configured evaluation defaults to the analysis start. Omit `[spinup]` for this default pathway: there are no artificial years or mandatory repeated-cycle gates. Incomplete calendar-year coverage is acceptable when the available daily prefix is complete from the chosen start. This is an initial-state assumption, not a claim of convergence from arbitrary profiles. See [the actual-2013 warm-up example](../examples/2013-warmup.toml). The repeated-year protocol below is an additional explicit option with its own complete-year and convergence requirements.
+For complete multi-year blocks, use the [calendar-block spin-up interface](multiyear-spinup.md).
+It keeps source forcing and output dates unexpanded. The `source_year` interface
+below remains available with its original artificial-prefix semantics.
+
+The usual Landsat 8/9 workflow runs alone. `forcing.peer_runs` is optional and should be set only for an explicitly requested comparison. For runs beginning in 2013, the default calibration scoring start is 2014-01-01, so available 2013 forcing supplies initialisation. Both calibration and normal runs (including `calibration.enabled = false`) start at physical field capacity, run the available 2013 record once and carry its state into 2014 without resetting. Resolved plans record the initialization and analysis windows; final SWB products preserve the actual calendar and mark warm-up days with `is_warmup`. Configured evaluation defaults to the analysis start. Omit `[spinup]` or leave `spinup.enabled = false` for this default pathway: there are no artificial years or mandatory repeated-cycle gates. Incomplete calendar-year coverage is acceptable when the available daily prefix is complete from the chosen start. This is an initial-state assumption, not a claim of convergence from arbitrary profiles. See [the actual-2013 warm-up example](../examples/2013-warmup.toml). The repeated-year protocol below is an additional explicit option with its own complete-year and convergence requirements.
 
 `forcing.peer_runs` lists the other run directories for a paired comparison. **Preprocess both runs before calibrating either.** Each run writes its own complete development forcing/physical-soil mask on both preparation grids. Initialization intersects corresponding masks with exact coordinates and CRS, preserving missing holdout forcing without using it to select cells. OpenLandMap's optional `soil.maximum_organic_matter_percent` also rejects hydraulic estimates outside the configured pedotransfer applicability domain.
 
 ```toml
 [spinup]
+enabled = true
 source_year = 2014
 repeats = 18
 tolerance = 0.005
diff --git a/docs/migration.md b/docs/migration.md
index 04b5e58..980ee34 100644
--- a/docs/migration.md
+++ b/docs/migration.md
@@ -1,5 +1,29 @@
 # Migration to the maintained package
 
+## Explicit repeated-spin-up switch (21 October 2026)
+
+Repeated forcing is disabled by default. Set `[spinup] enabled = true` to enable
+the existing single-year or multi-year protocol. An omitted or false `enabled`
+key leaves repetition disabled even when source years and a repetition count
+are present. This allows the settings to remain in a study configuration while
+the experiment is switched off. Actual run dates, physical field-capacity
+initialisation and explicit calibration windows remain unchanged.
+
+Update earlier repetition configurations by adding `enabled = true` when using
+this branch. Direct column calls likewise use
+`spinup={"enabled": True, "source_years": [2014, 2015], "repeats": 10}`.
+The existing component CLI/API `spinup_source_years` option is itself an explicit
+request to enable repetition. See [usage and disabled behaviour](multiyear-spinup.md).
+The daily physical equations and formulation identity are unchanged by this edit.
+
+## Depth-dependent extraction thresholds (21 October 2026)
+
+`swb-daily-depth-extraction-2p-v3` applies the frozen GSSM-to-WP dry ratio to extraction in 0–5 and 5–15 cm, uses the midpoint of adjusted and hydraulic WP in 15–30 cm, and retains WP below 30 cm. The continuous three-day supply relaxation, physical field-capacity initialisation, conservative flux update and zero-to-porosity physical bounds are retained. See the [equations, upper-bound assessment and limitations](swb-formulation.md).
+
+Use a new output directory and refit diffusivity and root beta. V2 fitted coefficients are rejected because the extraction equation changed. Existing immutable GSSM p01/p99 artifacts remain accepted without rewriting their hashes. Output diagnostics now include `below_reference_days_layer_*` for every layer, and calibration/product identities include `extraction_threshold_profile` and `physical_upper_bound`. The reference remains an empirical surface constraint; its depth transfer does not establish physiological validity or improved predictive skill.
+
+The depth-dependent physics also applies when repeated spin-up is enabled; the switch controls initialization only.
+
 ## September 2026 baseline
 
 GitHub and Gadi shared commit `e96879a`; the local checkout was four commits behind. Gadi additionally contained notebook/PBS improvements. Those two files were committed as `3fb0f23`, then GitHub, local and Gadi were fast-forwarded to that baseline before this update. This preserves the latest development rather than replacing it with an older checkout.
diff --git a/docs/multiyear-spinup.md b/docs/multiyear-spinup.md
new file mode 100644
index 0000000..4bec3d2
--- /dev/null
+++ b/docs/multiyear-spinup.md
@@ -0,0 +1,96 @@
+# Repeated calendar-block initialization
+
+Repeated spin-up is disabled by default. Set `spinup.enabled = true` to use a
+complete block at the beginning of the forcing record for repeated
+initialisation. For example:
+
+```toml
+[run]
+start = "2014-01-01"
+end = "2025-12-31"
+holdout_start = "2021-01-01"
+state_timing = "start"
+
+[calibration]
+enabled = true
+start = "2016-01-01"
+end = "2020-12-31"
+
+[spinup]
+enabled = true
+source_years = [2014, 2015]
+repeats = 10
+tolerance = 0.005
+minimum_fraction = 0.99
+convergence_mode = "wet_dry"
+convergence_action = "require"
+```
+
+This example integrates the 730-day block ten times (twenty simulated years).
+The count includes the retained 2014–2015 pass. The solver first integrates nine
+extra passes, discarding their trajectories, then continues through the actual
+2014–2025 record without resetting the soil state. Output dates remain unchanged.
+All 2014–2015 outputs are initialization and are excluded from fitting and
+evaluation. Calibration observations are never repeated.
+
+To disable repetition while retaining the other settings for later use, set:
+
+```toml
+[spinup]
+enabled = false
+```
+
+An omitted table, an empty table or an omitted `enabled` key also disables
+repetition. `enabled` must be a boolean. Disabled settings do not require source
+years, complete calendar blocks or convergence diagnostics. The ordinary run
+starts at physical field capacity and processes its actual dates once. Explicit
+calibration dates and any actual-forcing prefix remain in effect; the switch
+does not remove ordinary initialisation or rewrite the requested date range.
+
+Configurations created before the explicit switch need `enabled = true` to
+reproduce their repeated-year experiment with this version. Preserve their
+original source snapshot when reproducing an archived accepted run.
+
+`source_years` must be consecutive integer years forming a complete initial daily
+forcing block. Leap years retain every day. Calibration must begin after the
+block. The native workflow requires a development boundary before the holdout.
+The existing single-year `source_year` interface retains its earlier semantics.
+Do not specify both interfaces.
+
+The reference and Numba solvers share this implementation. Their public
+`soil_water_balance_1d` functions accept
+`spinup={"enabled": True, "source_years": [2014, 2015], "repeats": 10}`.
+Passing `None`, an empty mapping or `{"enabled": False}` disables repetition;
+omitting `enabled` in a supplied mapping also leaves it disabled.
+Gridded run and calibration
+APIs accept `spinup_source_years` and `spinup_repeats`; the command-line flags are
+`--spinup-source-years 2014 2015 --spinup-repeats 10`. Direct calibration also
+requires `--spinup-days 730` to exclude the retained block.
+Supplying `--spinup-source-years` explicitly enables repetition in these
+component APIs; omit the two repetition flags for the default disabled path.
+`--spinup-days` alone excludes an actual forcing prefix from calibration and
+does not repeat any forcing.
+
+Only one column's current cycle is allocated during initialization. Neither
+forcing grids nor spatial warm-up outputs are expanded with repetition count.
+Runtime still increases with the number of simulated days.
+
+## Interpreting convergence
+
+`wet_dry` compares empty and saturated initial states at the scoring boundary.
+`field_capacity_cycles` compares field-capacity starts with N and 2N total passes;
+passing it does not establish independence from other initial states. Each test
+requires every layer to meet the tolerance in the requested fraction of admitted
+cells. Default and newly calibrated parameters are checked separately.
+
+The default `convergence_action = "require"` stops on a failed diagnostic.
+An explicitly labelled fixed-duration research experiment may use `"record"`.
+This retains a failed receipt, emits a warning and records the interpretation in
+product initialization metadata. It does **not** turn the diagnostic into a pass
+or establish readiness for regional inference. Numerical conservation and
+calibration reproduction remain mandatory in either mode.
+
+Repeated forcing can reduce sensitivity to assumed initial water storage. It
+does not correct bias in forcing, hydraulic properties, rooting assumptions or
+the observation operator, and a longer initialization need not converge for
+every parameter setting.
diff --git a/docs/swb-formulation.md b/docs/swb-formulation.md
index b5b2a7f..738cee1 100644
--- a/docs/swb-formulation.md
+++ b/docs/swb-formulation.md
@@ -1,6 +1,6 @@
 # Conservative two-parameter SWB
 
-Formulation: `swb-daily-conservative-2p-v2`. Reference schema: `gssm-native-p01-p99-v1`.
+Formulation: `swb-daily-depth-extraction-2p-v3`. Reference schema: `gssm-native-p01-p99-v1`.
 
 The revised water balance separates physical storage, empirical surface-moisture references and root-extraction limits. Only the diffusivity scale `diff_factor` (10–10,000 mm) and the beta root-distribution coefficient `root_beta` (0.9–0.995) are calibrated. Root beta defaults to **0.96** in normal runs, the low-level root-distribution helper and the initial calibration seed; explicit values override this default. Physical saturation water content θs, wilting point WP, field capacity FC, saturated conductivity Ksat and the hydraulic exponent b remain fixed. The model layers end at 50, 150, 300, 600 and 1,000 mm. Layer 1 represents 0–5 cm and directly matches the GSSM observation depth.
 
@@ -10,7 +10,17 @@ For each admitted native GSSM cell, temporal p01 and p99 are calculated over 201
 
 Layer-1 hydraulic fields are calculated at the native soil-predictor cells before spatial averaging. Exact spherical overlap areas integrate these fields over the same union of complete GSSM cells. The implementation rejects incomplete native hydraulic coverage and invalid hydraulic values. This reference calculation preserves the nonlinear order of the pedotransfer calculation. Existing coarse and fine simulation soils retain their established preparation pathway; their spatial aggregation approximation remains a transfer limitation.
 
-Let Gd and Gw denote the area means of the native-cell p01 and p99. The fixed scalars are `c_d = min(1, Gd / mean(WP_1))` and `c_w = max(1, Gw / mean(θs_1))`. Shallow references are `d_i = c_d WP_i` and `w_i = c_w θs_i` for the first two layers. Deeper layers retain their unadjusted hydraulic references. Direct evaporation draws from layer 1; layer 2 supplies water through the hydraulic interface. The wet reference supports diagnostics and leaves physical storage capacity unchanged. A surface-derived scalar applied to layer 2 is an explicit shallow extrapolation.
+Let Gd and Gw denote the area means of the native-cell p01 and p99. The fixed scalars are `c_d = min(1, Gd / mean(WP_1))` and `c_w = max(1, Gw / mean(θs_1))`. The dry ratio is GSSM divided by hydraulic WP, so a lower GSSM reference reduces the extraction threshold. The first-layer reference is transferred to depth using the following fixed profile:
+
+| Layer | Depth | Extraction threshold d_i |
+|---|---|---|
+| 1 | 0–5 cm | `c_d WP_1` |
+| 2 | 5–15 cm | `c_d WP_2` |
+| 3 | 15–30 cm | `(WP_3 + c_d WP_3) / 2` |
+| 4 | 30–60 cm | `WP_4` |
+| 5 | 60–100 cm | `WP_5` |
+
+The third-layer midpoint combines its own hydraulic WP and its own fully adjusted value. For other layer grids, the adjustment weight is one over 0–150 mm, one-half over 150–300 mm and zero below 300 mm, averaged by layer thickness. Surface-derived thresholds in layers 2 and 3 are empirical extrapolations; GSSM provides no direct observation at those depths. Direct evaporation still draws from layer 1. The wet references remain `w_i = c_w θs_i` in the first two standard layers and support diagnostics only. They leave physical storage capacity unchanged.
 
 The JSON artifact is hashed over its complete canonical content. Set `[model].reference_file` in a configured workflow, or pass `reference_file` to calibration and simulation. Default and calibrated runs, initialisation checks and coarse-to-fine transfer use the same artifact. Native configured preparation retains the artifact's frozen complete-cell population and rejects loss of a member through invalid forcing or soil. Generic component runs without an artifact use the explicitly identified `physical-unadjusted` reference, with dry/wet factors of one.
 
@@ -18,7 +28,15 @@ The JSON artifact is hashed over its complete canonical content. Set `[model].re
 
 Physical state bounds are `0 ≤ θ_i ≤ θs_i`. Initialisation uses `FC = WP + available_water_capacity`. An initially dry layer receives no floor-restoration water. Root fractions integrate the beta profile over the actual rooted portion of each layer, using depth in centimetres and the fixed maximum root depth (1,000 mm by default). Fractions sum to one over the rooted profile, with no compensatory uptake.
 
-For the daily update, h = 1 day and `r = -expm1(-h/τ)`, with fixed `τ = 3 days`. The surface evaporation supply cap is `dz_1 θ_1²/(θ_1+d_1) r` in mm. Actual evaporation is the minimum of that cap and requested surface evaporation multiplied by h. The cap is zero at zero water and has a continuous limit when the dry reference is zero. Root extraction in layer i is bounded by both `root_fraction_i × requested_transpiration × h` and `dz_i max(θ_i−WP_i, 0) r`. WP constrains root extraction; evaporation and hydraulic redistribution can continue below WP and the dry reference. The closure is an empirical supply cap. Physiological interpretation requires additional evidence.
+For the daily update, h = 1 day and `r = -expm1(-h/τ)`, with fixed `τ = 3 days`. The surface evaporation supply cap is `dz_1 θ_1²/(θ_1+d_1) r` in mm. Actual evaporation is the minimum of that cap and requested surface evaporation multiplied by h. The cap is zero at zero water and has a continuous limit when the dry reference is zero. Root extraction in layer i is bounded by both `root_fraction_i × requested_transpiration × h` and `dz_i max(θ_i−d_i, 0) r`. This preserves the continuous three-day relaxation accepted on 18 October. With ample transpiration demand and no other fluxes, `θ_next = d_i + (θ_i − d_i) exp(−h/τ)` for states above d_i. Extraction stops at d_i; evaporation and hydraulic redistribution can continue below it. Initially lower states receive no restoration water. The older 1.25-times trigger and post-update floor restoration remain retired.
+
+Physical WP remains available for diagnostics and calculation of `FC = WP + available_water_capacity`. V3 deliberately replaces the WP extraction cutoff in the upper three layers with the empirical d_i profile. This changes the root-extraction assumption. A surface moisture percentile does not independently establish plant access to water below the hydraulic WP estimate; physiological interpretation and improved field skill require separate evidence.
+
+## Assessment of the upper and lower limits
+
+Hydraulic saturation/porosity remains the physical upper bound in every layer. Saturation represents filled pore space, as described in the [FAO drainage guidance](https://www.fao.org/4/R4082E/r4082e07.htm). A GSSM wet percentile is an observed climatological statistic. It can differ from saturation through sampling, scale and retrieval effects; it does not provide a direct estimate of pore volume. Therefore no GSSM multiplier expands or reduces the physical upper bound. Saturation excess continues through the conservative overflow pathway.
+
+WP is a plant-water extraction threshold and cannot represent a universal lower bound on total soil moisture. The [FAO soil-water discussion](https://www.fao.org/4/a0072e/a0072e07.htm) explicitly describes continued evaporation below permanent wilting point. The [FAO transpiration treatment](https://www.fao.org/4/x0490e/x0490e0e.htm) separately associates WP with cessation of root uptake. This distinction supports allowing soil moisture below WP, while the chosen reduction of the model extraction threshold remains an empirical development assumption. Physical model bounds remain zero to porosity; zero is a non-negativity limit and is not presented as a calibrated air-dry water content.
 
 Hydraulic saturation ratio is `S = θ/θs`, including the zero-water limit. Conductivity is `K = Ksat S^(2b+3)` in mm day⁻¹ and diffusivity is `D = diff_factor Ksat S^(b+2)` in mm² day⁻¹. The inherited interface discretisation uses the upper-layer coefficients: downward flux is `K_i + D_i(θ_i−θ_(i+1))/distance`. Bottom drainage is the smaller of the fixed drainage-slope times bottom-layer K and the fixed drainage upper limit. The lower drainage limit must be zero. An exactly empty upper layer has zero interface coefficients under this inherited discretisation; upward replenishment in that limiting state requires subsequent wetting. This is a structural limitation to retain in scientific interpretation.
 
@@ -30,11 +48,11 @@ This is a discrete daily water-balance formulation. Donor and receiver limits pr
 
 ## Outputs and calibration identity
 
-Every interval retains requested and actual E/T, unmet E/T, bottom drainage, saturation overflow, profile and layer mass residuals, signed interface transfers, layer transpiration and update counts (the compatibility variable `substeps` equals one on each valid day). Exposure diagnostics count days whose beginning-of-day state lies below WP or outside the shallow d/w reference envelope; they do not resolve within-day crossings. The empirical Gd/Gw envelope is distinct from the transferred d/w references and can be evaluated using the frozen artifact. The compatibility field `storage_adjustment` is identically zero.
+Every interval retains requested and actual E/T, unmet E/T, bottom drainage, saturation overflow, profile and layer mass residuals, signed interface transfers, layer transpiration and update counts (the compatibility variable `substeps` equals one on each valid day). Exposure diagnostics count beginning-of-day states below physical WP and below the extraction threshold in every layer, and above the wet reference in the first two layers; they do not resolve within-day crossings. Column metadata retains the full extraction-threshold and physical-WP arrays. Calibration and product identities record the extraction profile and physical upper-bound policy. The empirical Gd/Gw envelope is distinct from the transferred d/w references and can be evaluated using the frozen artifact. The compatibility field `storage_adjustment` is identically zero.
 
 The profile identity is `ΔW = P_eff − E_actual − T_actual − drainage − overflow`. Layer identities use the same interface volumes on either side. Missing forcing holds the state, sets the forcing-validity flag to zero and leaves interval budgets missing. Start/end state timestamps and the final restart state are retained.
 
-Calibration preserves the area-weighted native GSSM cell–day RMSE, fixed support, initialisation exclusion and coarse-to-fine transfer semantics. Seeds, optimizer populations, polishing, CSV output and downstream loading use exactly two fitted coefficients. Calibration CSVs and cache identities include formulation version, reference identity and fixed numerical settings. Subcycled-v1, four-parameter and unversioned calibrations raise explicit migration errors. Unmodified p01/p99 reference artifacts from v1 remain valid because their observations, hydraulic support and scalar definitions are independent of temporal integration; their hashes are preserved. Migration requires freezing the reference and recalibrating with this formulation; previous fits cannot be relabelled.
+Calibration preserves the area-weighted native GSSM cell–day RMSE, fixed support, initialisation exclusion and coarse-to-fine transfer semantics. Seeds, optimizer populations, polishing, CSV output and downstream loading use exactly two fitted coefficients. Calibration CSVs and cache identities include formulation version, reference identity and fixed numerical settings. Daily-v2, subcycled-v1, four-parameter and unversioned calibrations raise explicit migration errors. Unmodified p01/p99 reference artifacts from v1 and v2 remain valid because their observations, hydraulic support and scalar definitions are unchanged; their hashes are preserved. Migration requires recalibrating with this formulation and generating new trajectories; previous fits cannot be relabelled.
 
 ## Validation and interpretation
 
@@ -45,6 +63,6 @@ GSSM supplies reference statistics and calibration targets. In-situ observations
 
 ## Default initialisation with actual forcing
 
-Start each layer at physical field capacity on the first available simulation date. For a Landsat 8/9 record starting in 2013, in both calibration and normal runs, use the available actual 2013 daily forcing as warm-up, carry the resulting state into 2014 without resetting, and score from 2014-01-01. Do not prepend synthetic years. Omit the optional `[spinup]` table to select this pathway; its repeated-year mechanism is opt-in. The warm-up is excluded from the GSSM objective, while its forcing still affects subsequent states. Coverage must be checked. This protocol retains field capacity as an initial-state assumption and does not claim independence from every possible initial profile.
+Start each layer at physical field capacity on the first available simulation date. For a Landsat 8/9 record starting in 2013, in both calibration and normal runs, use the available actual 2013 daily forcing as warm-up, carry the resulting state into 2014 without resetting, and score from 2014-01-01. Do not prepend synthetic years. Omit the optional `[spinup]` table or leave `spinup.enabled = false` to select this pathway. Repetition requires `spinup.enabled = true`. The warm-up is excluded from the GSSM objective, while its forcing still affects subsequent states. Coverage must be checked. This protocol retains field capacity as an initial-state assumption and does not claim independence from every possible initial profile.
 
 Resolved plans and paired outputs record the initialization policy. The output `is_warmup` time flag separates the retained actual 2013 record from analysis dates, including when calibration is disabled. A normal run ending in 2013 is explicitly warm-up only. Later-starting studies and explicitly selected historical missions retain their requested dates.
diff --git a/examples/workflow.toml b/examples/workflow.toml
index 07d6f01..148fa70 100644
--- a/examples/workflow.toml
+++ b/examples/workflow.toml
@@ -9,7 +9,7 @@ center = [147.25, -35.05] # Default 5100 m square: 170 x 170 native 30 m pixels.
 # For a historical run, select its first supported date in 2013 or later.
 # If start is in 2013, calibration starts on 2014-01-01 unless explicitly set.
 # One daily SWB update; physical FC start; available 2013 is warm-up only.
-# No repeated years unless a [spinup] table is explicitly supplied.
+# Repeated spin-up is disabled by default; enable it with [spinup] enabled = true.
 grid_mode = "native"
 sensors = ["LC08", "LC09"] # Add "LT05", "LE07" explicitly when needed.
 processing_buffer_m = 6000
diff --git a/pysweb/swb/accelerated.py b/pysweb/swb/accelerated.py
index cb225a6..e9c3861 100644
--- a/pysweb/swb/accelerated.py
+++ b/pysweb/swb/accelerated.py
@@ -4,7 +4,7 @@
 Objective: Provide an optional numerically equivalent JIT backend for the default column solver.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-10-16
-Last updated: 2026-10-18
+Last updated: 2026-10-20
 Inputs: Daily forcing and column properties accepted by the maintained solver.
 Outputs: Identical state/budget schema, without silent reference fallback.
 Usage: from pysweb.swb.accelerated import soil_water_balance_1d
@@ -19,10 +19,10 @@
 
 def soil_water_balance_1d(precip_data, et_data, soil_properties, time, time_step=1.,
                           initial_soil_moisture=None, evap_fraction=None, diff_factor=1e3,
-                          transpiration_data=None, state_timing="start"):
+                          transpiration_data=None, state_timing="start", spinup=None):
     """Execute the compiled conservative kernel after common input validation."""
     result = _solve(precip_data, et_data, soil_properties, time, time_step,
         initial_soil_moisture, evap_fraction, diff_factor, transpiration_data,
-        state_timing, kernel=_kernel)
+        state_timing, kernel=_kernel, spinup=spinup)
     result.attrs["solver_backend"] = "numba"
     return result
diff --git a/pysweb/swb/calibrate.py b/pysweb/swb/calibrate.py
index e2b4f1e..6a93ccb 100644
--- a/pysweb/swb/calibrate.py
+++ b/pysweb/swb/calibrate.py
@@ -4,7 +4,7 @@
 Objective: Calibrate domain SWB parameters against reference surface soil moisture.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-04-19
-Last updated: 2026-10-20
+Last updated: 2026-10-21
 Inputs: Forcing and soil NetCDFs, reference SSM, calibration dates, bounds and random seed.
 Outputs: A CSV of converged parameters and fit diagnostics, with a provenance manifest.
 Usage: pysweb swb-calibrate --help
@@ -307,6 +307,8 @@ def _compute_rmse(
                     diff_factor=diff_factor,
                     transpiration_data=transpiration_vals[:, lat_idx, lon_idx],
                     state_timing=state_timing,
+                    **({"spinup": observation_context["spinup"]}
+                       if "spinup" in observation_context else {}),
                 )
             except (ValueError, FloatingPointError, ArithmeticError):
                 # Penalise this candidate instead of changing its calibration population.
@@ -431,6 +433,8 @@ def build_parser() -> argparse.ArgumentParser:
                         help="Area-weighted cell-day RMSE (default), or legacy RMSE of daily domain means.")
     parser.add_argument("--state-timing", choices=("start", "end"), default="start", help="Model state convention matched to reference observations.")
     parser.add_argument("--spinup-days", type=int, default=0, help="Initial simulated days excluded from the objective.")
+    parser.add_argument("--spinup-source-years", nargs="+", type=int, help="Complete initial years repeated before scoring.")
+    parser.add_argument("--spinup-repeats", type=int, default=1, help="Total passes including the retained source block.")
     parser.add_argument("--allow-unconverged", action="store_true", help="Explicitly retain a non-converged candidate with its status recorded.")
     parser.add_argument("--max-iter", type=int, default=30, help="Maximum iterations for optimization.")
     parser.add_argument("--surface-depth", type=float, default=50.0, help="Depth of surface observation (mm).")
@@ -784,6 +788,15 @@ def calibrate_domain(*, benchmark_output=None, polish_file=None, **kwargs):
         effective_precip_vals, et_vals, transpiration_vals, reference_vals,
         soil_valid, cell_weights, args.spinup_days, args.solver_backend, args.objective)
     observation_context["fixed_properties"] = {**reference_properties(args.reference_file), "max_root_depth": args.max_root_depth}
+    if args.spinup_source_years is not None:
+        from pysweb.swb.spinup import source_block
+        spinup = dict(enabled=True, source_years=args.spinup_source_years, repeats=args.spinup_repeats)
+        count, _ = source_block(time_index, spinup)
+        if args.spinup_days < count:
+            raise ValueError("Calibration must exclude the entire retained spin-up source block")
+        observation_context["spinup"] = spinup
+    elif args.spinup_repeats != 1:
+        raise ValueError("spinup_repeats requires spinup_source_years")
     if args.reference_file:
         from pysweb.swb.references import frozen_mask
         from pysweb.swb.formulation import load_reference
diff --git a/pysweb/swb/conservative.py b/pysweb/swb/conservative.py
index c74d4e3..441c862 100644
--- a/pysweb/swb/conservative.py
+++ b/pysweb/swb/conservative.py
@@ -3,7 +3,7 @@
 Objective: Advance a soil column with jointly limited conservative hydraulic fluxes and separate ET sinks.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-10-18
-Last updated: 2026-10-18
+Last updated: 2026-10-21
 Inputs: Validated float64 column properties, forcing and fixed closure settings.
 Outputs: States, interval budgets, interface fluxes, exposure durations and final state.
 Usage: Called by solver.py directly or compiled by accelerated.py.
@@ -25,7 +25,7 @@ def column_kernel(rain, et, tr, theta, por, wp, bcoef, ks, thick, dry, wet,
     """
     nt, nl = len(rain), len(theta)
     ns = min(2, nl)
-    count = 14 + nl - 1 + 3 * nl + 2 * ns
+    count = 14 + nl - 1 + 4 * nl + ns
     states = np.zeros((nt, nl))
     budgets = np.full((nt, count), np.nan)
     valid = np.zeros(nt, dtype=np.bool_)
@@ -43,7 +43,7 @@ def column_kernel(rain, et, tr, theta, por, wp, bcoef, ks, thick, dry, wet,
         qsum = np.zeros(max(0, nl-1))
         tsum = np.zeros(nl)
         below_wp = np.zeros(nl)
-        below_dry = np.zeros(ns)
+        below_dry = np.zeros(nl)
         above_wet = np.zeros(ns)
         esum, drainage, overflow = 0., 0., 0.
         h = dt
@@ -62,7 +62,9 @@ def column_kernel(rain, et, tr, theta, por, wp, bcoef, ks, thick, dry, wet,
         e = 0.
         if theta[0] > 0.:
             e = min(evap * h, thick[0] * theta[0] * theta[0] / (theta[0] + dry[0]) * release)
-        uptake = np.minimum(trans * roots * h, thick * np.maximum(theta - wp, 0.) * release)
+        # The continuous three-day supply relaxation now approaches the
+        # depth-dependent extraction threshold. WP remains a diagnostic.
+        uptake = np.minimum(trans * roots * h, thick * np.maximum(theta - dry, 0.) * release)
         drain = min(slope * cond[-1], dhigh) * h
         q = flux * h
         outgoing = uptake.copy()
@@ -97,9 +99,9 @@ def column_kernel(rain, et, tr, theta, por, wp, bcoef, ks, thick, dry, wet,
         for j in range(nl):
             if theta[j] < wp[j]:
                 below_wp[j] += h
-        for j in range(ns):
             if theta[j] < dry[j]:
                 below_dry[j] += h
+        for j in range(ns):
             if theta[j] > wet[j]:
                 above_wet[j] += h
         water -= uptake
@@ -146,6 +148,6 @@ def column_kernel(rain, et, tr, theta, por, wp, bcoef, ks, thick, dry, wet,
         offset += nl
         budgets[t, offset:offset+nl] = below_wp
         offset += nl
-        budgets[t, offset:offset+ns] = below_dry
-        budgets[t, offset+ns:offset+2*ns] = above_wet
+        budgets[t, offset:offset+nl] = below_dry
+        budgets[t, offset+nl:offset+nl+ns] = above_wet
     return states, budgets, valid, theta
diff --git a/pysweb/swb/formulation.py b/pysweb/swb/formulation.py
index a0a37c1..3885b64 100644
--- a/pysweb/swb/formulation.py
+++ b/pysweb/swb/formulation.py
@@ -3,7 +3,7 @@
 Objective: Define the versioned two-parameter SWB contract and reject incompatible calibrations.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-10-18
-Last updated: 2026-10-20
+Last updated: 2026-10-21
 Inputs: Fixed reference artifacts, model settings and calibration CSV files.
 Outputs: Validated fixed settings, artifact identities and fitted coefficients.
 Usage: import pysweb.swb.formulation
@@ -16,10 +16,12 @@
 import numpy as np
 import pandas as pd
 
-FORMULATION_VERSION = "swb-daily-conservative-2p-v2"
+FORMULATION_VERSION = "swb-daily-depth-extraction-2p-v3"
 TIME_INTEGRATION = "single-bounded-update"
+EXTRACTION_PROFILE = "gssm-wp-depth-blend-0-150-300mm-v1"
 # These observation-derived scalars are independent of the integration scheme.
-REFERENCE_COMPATIBLE_FORMULATIONS = {"swb-conservative-2p-v1", FORMULATION_VERSION}
+REFERENCE_COMPATIBLE_FORMULATIONS = {
+    "swb-conservative-2p-v1", "swb-daily-conservative-2p-v2", FORMULATION_VERSION}
 REFERENCE_VERSION = "gssm-native-p01-p99-v1"
 PARAMETER_NAMES = ("diff_factor", "root_beta")
 DEFAULT_ROOT_BETA = 0.96
@@ -66,6 +68,8 @@ def settings_identity(reference_file=None, *, max_root_depth=1000., drainage_slo
                       drainage_upper_limit=25., drainage_lower_limit=0.):
     reference = load_reference(reference_file)
     return dict(formulation_version=FORMULATION_VERSION, tau_days=TAU_DAYS,
+                extraction_threshold_profile=EXTRACTION_PROFILE,
+                physical_upper_bound="hydraulic_porosity",
                 time_integration=TIME_INTEGRATION, time_step_days=1.,
                 reference_artifact_id=reference["artifact_id"],
                 max_root_depth_mm=float(max_root_depth), drainage_slope=float(drainage_slope),
@@ -77,6 +81,27 @@ def reference_properties(reference_file=None):
     return {k: float(record[k]) for k in ("dry_reference_factor", "wet_reference_factor")}
 
 
+def extraction_thresholds(wilting_point, bottoms_mm, dry_factor):
+    """Transfer the surface dry ratio by depth, without changing hydraulic WP.
+
+    Full adjustment applies through 150 mm, half adjustment from 150 to 300 mm,
+    and none below 300 mm. Layers crossing these boundaries use thickness-
+    weighted adjustment, preserving the same depth rule on other layer grids.
+    """
+    wp = np.asarray(wilting_point, dtype=float)
+    bottoms = np.asarray(bottoms_mm, dtype=float)
+    if (wp.ndim != 1 or not len(wp) or wp.shape != bottoms.shape
+            or not np.isfinite(wp).all() or np.any(wp < 0)
+            or not np.isfinite(bottoms).all() or np.any(np.diff(np.r_[0., bottoms]) <= 0)
+            or not np.isfinite(dry_factor) or not 0 <= dry_factor <= 1):
+        raise ValueError("Invalid depth-dependent extraction threshold inputs")
+    tops = np.r_[0., bottoms[:-1]]
+    full = np.maximum(0., np.minimum(bottoms, 150.) - tops)
+    half = np.maximum(0., np.minimum(bottoms, 300.) - np.maximum(tops, 150.))
+    weight = (full + .5 * half) / (bottoms - tops)
+    return wp * (1. - (1. - dry_factor) * weight)
+
+
 def load_calibration(path, reference_file=None, **fixed_settings):
     """Never reinterpret a four-parameter or unversioned fit as the revised model."""
     # The daily trajectory can distinguish neighbouring float64 coefficients.
@@ -110,5 +135,5 @@ def budget_names(layers):
             + tuple(f"transpiration_layer_{i+1}" for i in range(layers))
             + tuple(f"layer_residual_{i+1}" for i in range(layers))
             + tuple(f"below_wp_days_layer_{i+1}" for i in range(layers))
-            + tuple(f"below_reference_days_layer_{i+1}" for i in range(min(2, layers)))
+            + tuple(f"below_reference_days_layer_{i+1}" for i in range(layers))
             + tuple(f"above_reference_days_layer_{i+1}" for i in range(min(2, layers))))
diff --git a/pysweb/swb/run.py b/pysweb/swb/run.py
index 5e90c8c..b8d9408 100644
--- a/pysweb/swb/run.py
+++ b/pysweb/swb/run.py
@@ -4,7 +4,7 @@
 Objective: Run serial or parallel gridded SWB simulations using the shared soil-water solver.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-04-17
-Last updated: 2026-10-20
+Last updated: 2026-10-21
 Inputs: Forcing and soil NetCDFs, simulation dates, model settings and optional calibration or parameter files.
 Outputs: A NetCDF of soil moisture, final states, water budgets and quality flags, plus a provenance manifest.
 Usage: pysweb swb-run --help
@@ -128,6 +128,8 @@ def build_parser() -> argparse.ArgumentParser:
     parser.add_argument("--drainage-lower-limit", type=float, default=0.0, help="Lower limit for drainage (mm day-1).")
     parser.add_argument("--solver-backend", choices=["reference", "numba"], default="reference", help="Optional verified accelerator for the standard solver equations.")
     parser.add_argument("--workers", type=int, default=1, help="Number of worker processes for spatial SWEB simulation.")
+    parser.add_argument("--spinup-source-years", nargs="+", type=int, help="Complete initial years repeated before scoring.")
+    parser.add_argument("--spinup-repeats", type=int, default=1, help="Total passes including the retained source block.")
     return parser
 
 
@@ -239,6 +241,7 @@ def _compute_row_soil_moisture(
     emit_warnings: bool = True,
     state_timing: str = "start",
     solver_backend: str = "reference",
+    spinup=None,
 ) -> np.ndarray:
     solver = soil_water_balance_1d
     if solver_backend == "numba":
@@ -293,6 +296,7 @@ def _compute_row_soil_moisture(
                 diff_factor = diff_factor,
                 transpiration_data = transpiration_series,
                 state_timing = state_timing,
+                **({"spinup": spinup} if spinup is not None else {}),
             )
         except Exception as exc:  # pragma: no cover - diagnostic output
             if emit_warnings:
@@ -324,6 +328,7 @@ def _init_process_run_worker(
     use_ndvi_root_depth: bool,
     state_timing: str,
     solver_backend: str = "reference",
+    spinup=None,
 ) -> None:
     global _PROCESS_RUN_STATE
     _PROCESS_RUN_STATE = {
@@ -339,6 +344,7 @@ def _init_process_run_worker(
         "use_ndvi_root_depth": bool(use_ndvi_root_depth),
         "state_timing": state_timing,
         "solver_backend": solver_backend,
+        "spinup": spinup,
     }
 
 
@@ -372,6 +378,7 @@ def _run_model_for_lat_process(lat_idx: int) -> tuple[int, np.ndarray]:
         emit_warnings = True,
         state_timing = state["state_timing"],
         solver_backend = state["solver_backend"],
+        spinup = state["spinup"],
     )
     return lat_idx, row_soil_moisture
 
@@ -390,6 +397,11 @@ def run_swb_workflow(**workflow_kwargs) -> Path | None:
         raise ValueError("Daily SWB forcing requires time_step=1 day.")
     if args.workers < 1:
         raise ValueError("workers must be >= 1")
+    spinup = None
+    if args.spinup_source_years is not None:
+        spinup = dict(enabled=True, source_years=args.spinup_source_years, repeats=args.spinup_repeats)
+    elif args.spinup_repeats != 1:
+        raise ValueError("spinup_repeats requires spinup_source_years")
 
     out_dir = Path(args.output_dir).expanduser().resolve()
     out_dir.mkdir(parents=True, exist_ok=True)
@@ -432,6 +444,9 @@ def run_swb_workflow(**workflow_kwargs) -> Path | None:
     except ValueError as exc:
         raise ValueError("Forcing grids must match exactly; preprocess onto one grid first.") from exc
     precip, effective_precip, et, transpiration = aligned[:4]
+    if spinup is not None:
+        from pysweb.swb.spinup import source_block
+        source_block(effective_precip.time.values, spinup)
     ndvi = aligned[4] if ndvi is not None else None
 
     precip = precip.transpose("time", args.lat_dim, args.lon_dim)
@@ -557,6 +572,7 @@ def run_swb_workflow(**workflow_kwargs) -> Path | None:
                 use_ndvi_root_depth = args.use_ndvi_root_depth,
                 state_timing = args.state_timing,
                 solver_backend = args.solver_backend,
+                spinup = spinup,
             )
             soil_moisture[:, :, lat_idx, :], final_soil_moisture[:, lat_idx, :], budget_values[:, :, lat_idx, :], forcing_valid[:, lat_idx, :] = row_soil_moisture
             processed_cells += n_lon
@@ -591,6 +607,7 @@ def run_swb_workflow(**workflow_kwargs) -> Path | None:
                 args.use_ndvi_root_depth,
                 args.state_timing,
                 args.solver_backend,
+                spinup,
             ),
         ) as executor:
             futures = {executor.submit(_run_model_for_lat_process, lat_idx) for lat_idx in range(n_lat)}
diff --git a/pysweb/swb/solver.py b/pysweb/swb/solver.py
index c2be798..700358e 100644
--- a/pysweb/swb/solver.py
+++ b/pysweb/swb/solver.py
@@ -3,7 +3,7 @@
 Objective: Solve daily soil-water balances and vertical redistribution for a layered soil column.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-04-17
-Last updated: 2026-10-20
+Last updated: 2026-10-21
 Inputs: Daily forcing, per-layer hydraulic properties, initial states and solver settings.
 Outputs: Layer moisture and final states, hydraulic coefficients, water-budget terms and quality flags.
 Usage: import pysweb.swb.solver
@@ -16,7 +16,8 @@
 
 from pysweb.contracts import layer_bottoms
 
-from pysweb.swb.formulation import (DEFAULT_ROOT_BETA, FORMULATION_VERSION, TAU_DAYS, TIME_INTEGRATION, budget_names, reject_legacy)
+from pysweb.swb.formulation import (DEFAULT_ROOT_BETA, EXTRACTION_PROFILE, FORMULATION_VERSION,
+    TAU_DAYS, TIME_INTEGRATION, budget_names, extraction_thresholds, reject_legacy)
 from pysweb.swb.conservative import column_kernel
 
 
@@ -39,7 +40,8 @@ def _resolve_sm_bounds(soil_properties):
 
 def _solve(precip_data, et_data, soil_properties, time, time_step=1.,
            initial_soil_moisture=None, evap_fraction=None, diff_factor=1e3,
-           transpiration_data=None, state_timing="start", *, kernel=column_kernel):
+           transpiration_data=None, state_timing="start", *, kernel=column_kernel,
+           spinup=None):
     """Validate once outside either numerical backend; retain missing-forcing states."""
     props = dict(soil_properties)
     reject_legacy(props)
@@ -86,9 +88,8 @@ def vector(value):
     wet_factor = float(props.get("wet_reference_factor", 1.))
     if not (np.isfinite(dry_factor) and 0 <= dry_factor <= 1 and np.isfinite(wet_factor) and wet_factor >= 1):
         raise ValueError("Invalid fixed dry/wet reference factors")
-    dry, wet = wp.copy(), por.copy()
+    dry, wet = extraction_thresholds(wp, bottoms, dry_factor), por.copy()
     shallow = bottoms <= 150.
-    dry[shallow] *= dry_factor
     wet[shallow] *= wet_factor
     roots = _compute_root_distribution_beta(props, beta=props.get("root_beta"),
         max_root_depth_mm=props.get("max_root_depth"), use_ndvi_root_depth=bool(props.get("use_ndvi_root_depth", False)))
@@ -100,6 +101,25 @@ def vector(value):
     tau = float(props.get("supply_tau_days", TAU_DAYS))
     if not (np.isfinite(tau) and tau > 0) or not np.isfinite(slope) or slope < 0 or not np.isfinite(upper) or upper < 0:
         raise ValueError("Invalid drainage or fixed closure settings")
+    original_initial = initial.copy()
+    spinup_metadata = None
+    from pysweb.swb.spinup import enabled_spinup
+    spinup = enabled_spinup(spinup)
+    if spinup is not None:
+        from pysweb.swb.spinup import source_block
+        if time_step != 1.:
+            raise ValueError("Repeated calendar spin-up requires daily forcing")
+        count, spinup_metadata = source_block(times, spinup)
+        for _ in range(spinup_metadata["repeats"] - 1):
+            _, _, warm_valid, initial = kernel(
+                rain[:count], et[:count], tr[:count], initial.copy(), por, wp, b, ks,
+                thick, dry, wet, roots, float(diff_factor), float(time_step),
+                slope, upper, tau, False)
+            if not np.all(warm_valid):
+                raise ValueError("Repeated spin-up forcing contains invalid values")
+        if not (np.isfinite(rain[:count]).all() and np.isfinite(et[:count]).all()
+                and np.isfinite(tr[:count]).all()):
+            raise ValueError("Repeated spin-up forcing contains invalid values")
     states, budgets, valid, final = kernel(rain, et, tr, initial.copy(), por, wp, b, ks,
         thick, dry, wet, roots, float(diff_factor), float(time_step), slope, upper,
         tau, state_timing == "end")
@@ -108,16 +128,21 @@ def vector(value):
         final_state=final, final_state_time=times[-1]+pd.Timedelta(days=time_step),
         formulation_version=FORMULATION_VERSION, supply_tau_days=tau, time_integration=TIME_INTEGRATION, time_step_days=float(time_step),
         water_balance={key: budgets[:, i] for i, key in enumerate(budget_names(len(bottoms)))},
-        forcing_valid=valid, dry_reference=dry, wet_reference=wet)
+        forcing_valid=valid, dry_reference=dry, wet_reference=wet,
+        extraction_threshold=dry.copy(), extraction_threshold_profile=EXTRACTION_PROFILE,
+        hydraulic_wilting_point=wp.copy(), physical_upper_bound="hydraulic_porosity")
+    if spinup_metadata is not None:
+        result.attrs.update(spinup=spinup_metadata, pre_spinup_initial_state=original_initial)
     return result
 
 
 def soil_water_balance_1d(precip_data, et_data, soil_properties, time, time_step=1.,
                           initial_soil_moisture=None, evap_fraction=None, diff_factor=1e3,
-                          transpiration_data=None, state_timing="start"):
+                          transpiration_data=None, state_timing="start", spinup=None):
     """One conservative update per interval; daily by default, initialised at FC."""
     return _solve(precip_data, et_data, soil_properties, time, time_step,
-        initial_soil_moisture, evap_fraction, diff_factor, transpiration_data, state_timing)
+        initial_soil_moisture, evap_fraction, diff_factor, transpiration_data,
+        state_timing, spinup=spinup)
 
 _DEFAULT_NDVI_ROOT_QUANTILES = (0.2, 0.4, 0.6, 0.8, 1.0)
 _DEFAULT_NDVI_ROOT_DEPTHS_CM = (5.0, 15.0, 30.0, 60.0, 100.0)
diff --git a/pysweb/swb/spinup.py b/pysweb/swb/spinup.py
index f5b410d..ce47ad4 100644
--- a/pysweb/swb/spinup.py
+++ b/pysweb/swb/spinup.py
@@ -3,7 +3,7 @@
 Objective: Assemble development-only forcing cycles and verify initial-state convergence.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-10-17
-Last updated: 2026-10-18
+Last updated: 2026-10-21
 Inputs: Daily forcing calendars, development boundary, hydraulic properties and parameters.
 Outputs: Explicit cycle mappings and wet/dry differences on fixed admitted cells.
 Usage: from pysweb.swb.spinup import cycle_indices, verify_convergence
@@ -14,6 +14,47 @@
 import xarray as xr
 
 
+def enabled_spinup(configuration):
+    """Return an explicitly enabled configuration; repetition is off by default."""
+    if configuration is None:
+        return None
+    if not isinstance(configuration, dict):
+        raise ValueError("spinup must be a configuration mapping or None")
+    enabled = configuration.get("enabled", False)
+    if not isinstance(enabled, bool):
+        raise ValueError("spinup.enabled must be a boolean")
+    return dict(configuration) if enabled else None
+
+
+def source_block(dates, configuration):
+    """Validate a complete initial calendar block; repetitions include its retained pass.
+
+    Only the first ``repeats - 1`` passes are discarded. The last pass remains
+    on its original dates and must be excluded from calibration/evaluation.
+    """
+    years = configuration.get("source_years")
+    repeats = configuration.get("repeats")
+    if "source_year" in configuration:
+        raise ValueError("Specify source_years or legacy source_year, not both")
+    if (not isinstance(years, (list, tuple)) or not years
+            or any(not isinstance(y, int) or isinstance(y, bool) for y in years)
+            or list(years) != list(range(years[0], years[-1] + 1))):
+        raise ValueError("source_years must be consecutive integer calendar years")
+    if not isinstance(repeats, int) or isinstance(repeats, bool) or repeats < 1:
+        raise ValueError("repeats must be a positive integer")
+    dates = pd.DatetimeIndex(dates)
+    expected = pd.date_range(f"{years[0]}-01-01", f"{years[-1]}-12-31")
+    if not dates[:len(expected)].equals(expected):
+        raise ValueError("Complete source years must form the initial daily forcing block")
+    return len(expected), dict(
+        source_years=list(years), repeats=repeats, cycle_days=len(expected),
+        initialization_days=len(expected) * repeats,
+        artificial_days=len(expected) * (repeats - 1),
+        scoring_start_minimum=f"{years[-1] + 1}-01-01",
+        repeat_count_includes_retained_source_block=True,
+        limitation="Repeated forcing initializes numerics; it does not reconstruct historical moisture")
+
+
 def cycle_indices(dates, actual_start, repeats, source_year, development_end):
     """Repeat a complete non-leap development year ending in the correct phase."""
     dates = pd.DatetimeIndex(dates)
@@ -96,7 +137,10 @@ def verify_convergence(forcing, soil, mask, parameters, *, scoring_start,
         raise ValueError("Complete daily initialization history required")
     if mode not in {"wet_dry", "field_capacity_cycles"}:
         raise ValueError("Unknown initialization convergence mode")
-    if mode == "field_capacity_cycles" and (
+    compact = enabled_spinup((model_options or {}).get("spinup"))
+    if compact:
+        _, compact_metadata = source_block(times, compact)
+    if mode == "field_capacity_cycles" and not compact and (
             not isinstance(cycle_days, int) or isinstance(cycle_days, bool)
             or cycle_days < 365 or cycle_days % 365 or cycle_days > len(times)):
         raise ValueError("Field-capacity comparison requires complete repeated-year cycle_days")
@@ -106,6 +150,7 @@ def verify_convergence(forcing, soil, mask, parameters, *, scoring_start,
     options = dict(drainage_slope=.5, drainage_upper_limit=25., drainage_lower_limit=0.,
                    use_ndvi_root_depth=False)
     options.update(model_options or {})
+    options["spinup"] = compact
     if options.get("nan_to_zero"):
         raise ValueError("Initialization cannot replace missing forcing with zero")
     values = {k: history[k].transpose("time", *spatial).values for k in ("effective_precip", "et", "E_t")}
@@ -151,6 +196,9 @@ def verify_convergence(forcing, soil, mask, parameters, *, scoring_start,
             max_root_depth=options.get("max_root_depth", 1000.),
             **{k: options[k] for k in ("drainage_slope", "drainage_upper_limit", "drainage_lower_limit")}),
         parameters=dict(parameters), failure_action="Extend development-only initialization and refit")
+    if compact:
+        record.update(repeated_source_block=compact_metadata,
+                      initialization_days=len(times)+compact_metadata["artificial_days"])
     data = xr.DataArray(difference, dims=("layer",)+spatial,
         coords={"layer": soil.layer, **{d: mask[d] for d in spatial}}, name="initialization_difference")
     data.attrs["units"] = "m3 m-3"
@@ -183,14 +231,17 @@ def _convergence_cell(cell, context=None):
     for index, initial in enumerate(initial_states):
         series = {key: values[key][:, y, x] for key in ("effective_precip", "et", "E_t")}
         dates = times
-        if mode == "field_capacity_cycles" and index == 1:
+        spinup = options.get("spinup")
+        if mode == "field_capacity_cycles" and index == 1 and spinup:
+            spinup = {**spinup, "repeats": 2 * spinup["repeats"]}
+        elif mode == "field_capacity_cycles" and index == 1:
             # Duplicate only the artificial prefix; append the real prefix once.
             series = {key: np.concatenate([value[:cycle_days], value]) for key, value in series.items()}
             dates = pd.date_range(end=times[-1], periods=len(times)+cycle_days)
         result = soil_water_balance_1d(series["effective_precip"],
             series["et"], props, dates, initial_soil_moisture=initial,
             diff_factor=parameters["diff_factor"], transpiration_data=series["E_t"],
-            state_timing="start")
+            state_timing="start", **({"spinup": spinup} if spinup else {}))
         if not np.all(result.attrs["forcing_valid"]):
             raise ValueError("Solver rejected initialization forcing")
         final.append(result.attrs["final_state"])
diff --git a/pysweb/workflow.py b/pysweb/workflow.py
index ee23196..1045d7e 100644
--- a/pysweb/workflow.py
+++ b/pysweb/workflow.py
@@ -3,7 +3,7 @@
 Objective: Run explicitly configured stages with shared paths and model settings.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-09-12
-Last updated: 2026-10-18
+Last updated: 2026-10-21
 Inputs: Study TOML configuration, requested stages and optional dry-run flag.
 Outputs: Resolved configuration and products from the selected model stages.
 Usage: pysweb workflow --config examples/workflow.toml --dry-run
@@ -185,7 +185,8 @@ def load_config(path):
             raise ValueError(f"Unknown {section} settings: {sorted(unknown)}")
     for section, keys in {
         "landsat": {"manifest_file", "sensor_windows", "wrs", "gapfill"},
-        "spinup": {"source_year", "repeats", "tolerance", "minimum_fraction", "solver_backend", "convergence_mode"},
+        "spinup": {"enabled", "source_year", "source_years", "repeats", "tolerance", "minimum_fraction",
+                   "solver_backend", "convergence_mode", "convergence_action"},
         "forcing": {"peer_runs", "reconstruction", "precipitation_qc", "retrieval"},
         "evaluation": {"longitude", "latitude", "depth_interval_mm", "observations_file", "regimes_file",
                        "regime_column", "pairs", "phases", "scoring_start"},
@@ -209,7 +210,9 @@ def load_config(path):
     for key, raw in retrieval.get('meteorology_files', {}).items():
         value = Path(raw).expanduser()
         retrieval['meteorology_files'][key] = str(value if value.is_absolute() else (path.parent/value).resolve())
-    if run['grid_mode'] == 'legacy' and any(config.get(k) for k in ('landsat','spinup','forcing','evaluation')):
+    from pysweb.swb.spinup import enabled_spinup
+    if run['grid_mode'] == 'legacy' and (enabled_spinup(config.get('spinup'))
+            or any(config.get(k) for k in ('landsat','forcing','evaluation'))):
         raise ValueError('Configured protocol integration requires grid_mode=native')
     if run['grid_mode'] == 'native':
         from pysweb.workflow_native import native_plan
diff --git a/pysweb/workflow_native.py b/pysweb/workflow_native.py
index dc5abaa..7c3afb0 100644
--- a/pysweb/workflow_native.py
+++ b/pysweb/workflow_native.py
@@ -3,7 +3,7 @@
 Objective: Orchestrate native downloads, GSSM calibration and shared Landsat-grid products.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-10-17
-Last updated: 2026-10-18
+Last updated: 2026-10-21
 Inputs: Validated workflow configuration and selected execution stages.
 Outputs: Two final NetCDFs, coarse calibration inputs, parameters and grid checks.
 Usage: pysweb workflow --config examples/workflow.toml
@@ -18,6 +18,7 @@
 from pysweb.io.grids import open_template, subset_template, centered_template
 from pysweb.io.products import crop_product, assert_identical_grids
 from pysweb.io.provenance import atomic_json
+from pysweb.swb.spinup import enabled_spinup
 
 
 def buffered_extent(extent, metres):
@@ -82,16 +83,30 @@ def native_plan(config, stages):
         raise ValueError('Native calibration requires the area-weighted cell-day objective')
     if config.get('model',{}).get('nan_to_zero'):
         raise ValueError('Native workflow preserves missing forcing; nan_to_zero would change calibration/run semantics')
-    spinup = config.get('spinup')
+    spinup = enabled_spinup(config.get('spinup'))
+    compact = bool(spinup and 'source_years' in spinup)
+    compact_metadata = None
     if spinup:
         from pysweb.swb.spinup import cycle_indices
         if not boundary:
             raise ValueError('Configured spin-up requires an explicit development boundary')
-        cycle_indices(pd.date_range(start, end), start, spinup['repeats'], spinup['source_year'], boundary)
+        if compact:
+            from pysweb.swb.spinup import source_block
+            _, compact_metadata = source_block(pd.date_range(start, end), spinup)
+            if not enabled and 'start' not in calibration:
+                first = pd.Timestamp(compact_metadata['scoring_start_minimum'])
+            if first < pd.Timestamp(compact_metadata['scoring_start_minimum']):
+                raise ValueError('Scoring must begin after the retained spin-up source block')
+            if pd.Timestamp(compact_metadata['scoring_start_minimum']) > boundary:
+                raise ValueError('Spin-up source years must precede the development boundary')
+        else:
+            cycle_indices(pd.date_range(start, end), start, spinup['repeats'], spinup['source_year'], boundary)
         if not 0 < spinup['tolerance'] or not 0 < spinup['minimum_fraction'] <= 1:
             raise ValueError('Invalid spin-up acceptance thresholds')
         if spinup.get('convergence_mode', 'wet_dry') not in {'wet_dry', 'field_capacity_cycles'}:
             raise ValueError('Invalid spin-up convergence_mode')
+        if spinup.get('convergence_action', 'require') not in {'require', 'record'}:
+            raise ValueError('Invalid spin-up convergence_action')
         if run['state_timing'] != 'start':
             raise ValueError('Verified spin-up workflow requires start-of-day states')
     landsat = config.get('landsat', {})
@@ -99,16 +114,26 @@ def native_plan(config, stages):
         raise ValueError('L7 gap-fill donors must end within development')
     if config.get('forcing') and not boundary:
         raise ValueError('Forcing support requires a development boundary')
-    analysis_start = first if enabled else default_scoring_start
+    analysis_start = first if enabled or compact else default_scoring_start
     warmup_end = min(end, analysis_start-pd.Timedelta(days=1))
     warmup = [str(start.date()), str(warmup_end.date())] if warmup_end >= start else None
     initialization = dict(initial_state='physical_field_capacity',
-        mode='repeated_year' if spinup else 'actual_forcing_prefix',
+        repeated_spinup_enabled=bool(spinup),
+        mode='repeated_source_block' if compact else ('repeated_year' if spinup else 'actual_forcing_prefix'),
         actual_start=str(start.date()), warmup_window=warmup,
         analysis_start=str(analysis_start.date()),
         actual_warmup_days=max(0, (warmup_end-start).days+1),
-        artificial_warmup_days=365*spinup['repeats'] if spinup else 0,
+        artificial_warmup_days=(compact_metadata['artificial_days'] if compact
+                               else (365*spinup['repeats'] if spinup else 0)),
         interpretation='Field-capacity assumption; actual prefix alone does not establish initial-state independence')
+    if compact:
+        initialization['repeated_source_block'] = compact_metadata
+    if spinup:
+        initialization['convergence_action'] = spinup.get('convergence_action', 'require')
+        initialization['interpretation'] = (
+            'Fixed-duration initialization experiment; inspect convergence receipts before interpreting results'
+            if spinup.get('convergence_action') == 'record' else
+            'Repeated forcing; convergence must pass the explicitly selected diagnostic')
     base = Path(run['root'])/run['name']
     return {'stages':list(stages),'run_directory':str(base), 'windows':[[run['start'],run['end']]],
             'calibration_scoring_window':[str(first.date()),str(last.date())] if enabled else None,
@@ -127,6 +152,7 @@ def run_native(config, stages, dry_run=False):
     from pysweb.swb.paired import run_default_and_calibrated
 
     plan = native_plan(config, stages)
+    spinup = enabled_spinup(config.get('spinup'))
     if dry_run:
         return plan
     run, soil = config['run'], config.get('soil', {'source':'openlandmap'})
@@ -273,7 +299,7 @@ def run_native(config, stages, dry_run=False):
                 reference_file=str(reference),reference_var=soil.get('reference_var','reference_ssm'))
         pysweb.swb.preprocess(**common,end_date=end,grid_mode='landsat',grid_template=str(seb_final),
             grid_template_var='ET',output_dir=str(fine),skip_reference_ssm=True)
-        if config.get('spinup') or config.get('forcing'):
+        if spinup or config.get('forcing'):
             from pysweb.workflow_protocol import write_support
             if enabled:
                 write_support(coarse, start, fit_end, boundary)
@@ -290,8 +316,7 @@ def forcing(directory,last,projected=False):
             lat_dim='y' if projected else 'lat',lon_dim='x' if projected else 'lon',
             state_timing=run['state_timing'],workers=run['workers'])
 
-    protocol = bool(config.get('spinup') or config.get('forcing'))
-    spinup = config.get('spinup')
+    protocol = bool(spinup or config.get('forcing'))
     def initialize(directory, last, projected=False):
         from pysweb.workflow_protocol import initialize_inputs, initialized_forcing
         peers = [Path(p)/directory.name for p in config.get('forcing', {}).get('peer_runs', [])]
@@ -373,11 +398,11 @@ def initialize(directory, last, projected=False):
                     values = {k:float(parameters[k]) for k in defaults}
                     convergence_gate(fine,record,mask,values,spinup,'calibrated',model,workers=run['workers'])
         output = run_default_and_calibrated(**args,**model, precip=rain,output_dir=str(products),
-            output_file='SWEB_with_spinup.nc' if spinup else 'SWEB_RZSM.nc',
+            output_file='SWEB_with_spinup.nc' if spinup and 'source_years' not in spinup else 'SWEB_RZSM.nc',
             calibration_file=calibration_file if enabled else None,
             initialization=plan['initialization'],
             auxiliary_inputs={'soil_evaporation_forcing':evaporation})
-        if spinup:
+        if spinup and 'source_years' not in spinup:
             from pysweb.workflow_protocol import trim_initialization
             output = trim_initialization(output, products/'SWEB_RZSM.nc', start)
         if protocol:
diff --git a/pysweb/workflow_protocol.py b/pysweb/workflow_protocol.py
index 348166b..ceb95b8 100644
--- a/pysweb/workflow_protocol.py
+++ b/pysweb/workflow_protocol.py
@@ -3,7 +3,7 @@
 Objective: Connect configured support, initialization and evaluation gates to native runs.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-10-17
-Last updated: 2026-10-18
+Last updated: 2026-10-21
 Inputs: Native prepared directories, paired workflow paths and scientific protocol settings.
 Outputs: Frozen support, explicit initialization files, convergence receipts and evaluation.
 Usage: Internal orchestration for pysweb.workflow_native.
@@ -63,7 +63,15 @@ def write_support(directory, start, end, development_end):
 def initialize_inputs(directory, start, end, scoring_start, development_end,
                       spinup=None, peer_directories=()):
     """Freeze support and retain actual forcing; prepend cycles only when requested."""
-    from pysweb.swb.spinup import prepend_cycles
+    from pysweb.swb.spinup import enabled_spinup, prepend_cycles
+    spinup = enabled_spinup(spinup)
+    compact = bool(spinup and "source_years" in spinup)
+    compact_metadata = None
+    if compact:
+        from pysweb.swb.spinup import source_block
+        _, compact_metadata = source_block(pd.date_range(start, end), spinup)
+        if pd.Timestamp(scoring_start) < pd.Timestamp(compact_metadata["scoring_start_minimum"]):
+            raise ValueError("Scoring must exclude the retained spin-up source block")
     directory = Path(directory)
     masks = []
     for folder in [directory, *map(Path, peer_directories)]:
@@ -84,12 +92,12 @@ def initialize_inputs(directory, start, end, scoring_start, development_end,
     actual_tag = f'{start.replace("-", "")}_{end.replace("-", "")}'
     simulation_start = start
     mapping = {}
-    if spinup:
+    if spinup and not compact:
         simulation_start = str((pd.Timestamp(start)-pd.Timedelta(days=365*spinup["repeats"])).date())
     for name, (stem, variable) in FORCING_FILES.items():
         source = directory/f"{stem}_daily_{actual_tag}.nc"
         destination = prepared/source.name
-        if spinup:
+        if spinup and not compact:
             prepend_cycles(source, destination, start, end, source_year=spinup["source_year"],
                 repeats=spinup["repeats"], development_end=development_end)
         else:
@@ -98,7 +106,7 @@ def initialize_inputs(directory, start, end, scoring_start, development_end,
     reference = directory/f"reference_ssm_daily_{actual_tag}.nc"
     if reference.exists():
         destination = prepared/reference.name
-        if spinup:
+        if spinup and not compact:
             prepend_cycles(reference, destination, start, end, source_year=spinup["source_year"],
                 repeats=spinup["repeats"], development_end=development_end,
                 reference=True, scoring_start=scoring_start)
@@ -109,8 +117,10 @@ def initialize_inputs(directory, start, end, scoring_start, development_end,
         scoring_start=scoring_start, development_end=development_end,
         spinup_days=(pd.Timestamp(scoring_start)-pd.Timestamp(simulation_start)).days,
         spinup=spinup, initial_state="physical_field_capacity",
-        initialization_mode="repeated_year" if spinup else "actual_forcing_prefix",
-        artificial_warmup_days=(pd.Timestamp(start)-pd.Timestamp(simulation_start)).days,
+        repeated_spinup_enabled=bool(spinup),
+        initialization_mode="repeated_source_block" if compact else ("repeated_year" if spinup else "actual_forcing_prefix"),
+        artificial_warmup_days=(compact_metadata["artificial_days"] if compact else
+                               (pd.Timestamp(start)-pd.Timestamp(simulation_start)).days),
         cells=int(mask.sum()), peers=list(map(str, peer_directories)), files=mapping)
     sources = [directory/f'{stem}_daily_{actual_tag}.nc' for stem, _ in FORCING_FILES.values()]
     sources += [directory/name for name in SOIL_FILE_STEMS.values()]
@@ -130,6 +140,9 @@ def initialized_forcing(record, directory, *, workers, state_timing, projected=F
         state_timing=state_timing, lat_dim="y" if projected else "lat", lon_dim="x" if projected else "lon")
     for name, argument in keys.items():
         result[argument], result[argument+"_var"] = record["files"][name]
+    if record.get("initialization_mode") == "repeated_source_block":
+        result.update(spinup_source_years=record["spinup"]["source_years"],
+                      spinup_repeats=record["spinup"]["repeats"])
     return result
 
 
@@ -147,16 +160,28 @@ def convergence_gate(directory, record, mask, parameters, config, label, model_o
     soil = xr.Dataset(arrays).assign_coords(layer_depth=("layer", infer_layer_bottoms(arrays, None)))
     # Dataset alignment expands the other histories to NDVI's calendar; verification
     # explicitly selects the initialization prefix again.
+    options = dict(model_options or {})
+    if record.get("initialization_mode") == "repeated_source_block":
+        options["spinup"] = record["spinup"]
     gate, differences = verify_convergence(xr.Dataset(fields), soil, mask, parameters,
         scoring_start=record["scoring_start"], development_end=record["development_end"],
         tolerance=config["tolerance"], minimum_fraction=config["minimum_fraction"],
-        solver_backend=config.get("solver_backend", "numba"), model_options=model_options, workers=workers,
+        solver_backend=config.get("solver_backend", "numba"), model_options=options, workers=workers,
         mode=config.get("convergence_mode", "wet_dry"),
         cycle_days=(pd.Timestamp(record["actual_start"])-pd.Timestamp(record["simulation_start"])).days)
+    action = config.get("convergence_action", "require")
+    if action not in {"require", "record"}:
+        raise ValueError("Invalid spin-up convergence_action")
+    gate["convergence_action"] = action
+    if action == "record":
+        gate["failure_action"] = "Retain fixed-duration experiment with failed convergence explicitly reported; no convergence endorsement"
     atomic_json(directory/"initialized"/f"{label}_convergence.json", gate)
     atomic_netcdf(differences, directory/"initialized"/f"{label}_convergence.nc")
-    if not gate["passed"]:
+    if not gate["passed"] and action == "require":
         raise ValueError(f"Initialization convergence failed ({gate['mode']}); extend development-only cycles and refit")
+    if not gate["passed"]:
+        import warnings
+        warnings.warn(f"{label} initialization did not converge; retaining explicitly configured fixed-duration experiment", RuntimeWarning)
     return gate
 
 
diff --git a/tests/swb/test_conservative_formulation.py b/tests/swb/test_conservative_formulation.py
index 9ec825d..52d88dd 100644
--- a/tests/swb/test_conservative_formulation.py
+++ b/tests/swb/test_conservative_formulation.py
@@ -3,7 +3,7 @@
 Objective: Independently verify physical storage, analytical limits, conservative budgets and migration.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-10-18
-Last updated: 2026-10-20
+Last updated: 2026-10-21
 Inputs: Synthetic columns, extreme forcing and fitted float64 coefficient fixtures.
 Outputs: Numerical assertions for the revised SWB contract.
 Usage: python -m pytest tests/swb/test_conservative_formulation.py
@@ -49,7 +49,7 @@ def test_evaporation_crosses_wp_and_reference_without_adding_water():
     assert np.all(np.diff(trace) < 0) and 0 < trace[-1] < .01
     assert trace[0] < .15 and np.any(trace < .06)
     assert np.all(np.diff(result.attrs["water_balance"]["actual_evaporation"]) < 0)
-    for initial in [0., .02, .06, .14]:
+    for initial in [0., .02, .06]:
         no_sink = soil_water_balance_1d([0.], [100.], p, ["2020-01-01"],
             transpiration_data=[100.], initial_soil_moisture=[initial])
         np.testing.assert_array_equal(no_sink.attrs["final_state"], [initial])
@@ -114,10 +114,11 @@ def test_subcycling_settings_and_prior_calibrations_are_rejected(tmp_path):
         soil_water_balance_1d([0.], [0.], props(max_substep_days=1/768), ["2020-01-01"])
     path = tmp_path / "subcycled.csv"
     identity = settings_identity()
-    identity["formulation_version"] = "swb-conservative-2p-v1"
-    pd.DataFrame([dict(diff_factor=1000., root_beta=.96, **identity)]).to_csv(path, index=False)
-    with pytest.raises(ValueError, match="formulation_version"):
-        load_calibration(path)
+    for previous in ["swb-conservative-2p-v1", "swb-daily-conservative-2p-v2"]:
+        identity["formulation_version"] = previous
+        pd.DataFrame([dict(diff_factor=1000., root_beta=.96, **identity)]).to_csv(path, index=False)
+        with pytest.raises(ValueError, match="formulation_version"):
+            load_calibration(path)
 
 
 @pytest.mark.parametrize("initial", [0., .01, .149, .449, .45])
diff --git a/tests/swb/test_extraction_thresholds.py b/tests/swb/test_extraction_thresholds.py
new file mode 100644
index 0000000..db7c787
--- /dev/null
+++ b/tests/swb/test_extraction_thresholds.py
@@ -0,0 +1,132 @@
+"""
+Script: test_extraction_thresholds.py
+Objective: Verify depth-dependent extraction, continuous relaxation and conserved water independently.
+Author: Yi Yu (with assistance from Codex)
+Created: 2026-10-21
+Last updated: 2026-10-21
+Inputs: Analytical drydowns and synthetic layered columns with known hydraulic properties.
+Outputs: Assertions for threshold transfer, extraction, physical bounds and both numerical backends.
+Usage: python -m pytest tests/swb/test_extraction_thresholds.py
+Dependencies: numpy, pandas, pytest; pysweb and optional numba.
+"""
+import numpy as np
+import pandas as pd
+import pytest
+
+
+@pytest.fixture(params=["reference", "numba"])
+def solve(request):
+    if request.param == "numba":
+        pytest.importorskip("numba")
+        from pysweb.swb.accelerated import soil_water_balance_1d
+    else:
+        from pysweb.swb.solver import soil_water_balance_1d
+    return soil_water_balance_1d
+
+
+def properties():
+    return dict(layer_depth=np.array([50., 150., 300., 600., 1000.]),
+        layer_thickness=np.array([50., 100., 150., 300., 400.]),
+        porosity=np.full(5, .5), wilting_point=np.array([.15, .18, .20, .25, .30]),
+        available_water_capacity=np.full(5, .05), b_coefficient=np.ones(5),
+        conductivity_sat=np.zeros(5), root_beta=.96, dry_reference_factor=.4)
+
+
+def test_all_layers_follow_analytical_three_day_extraction_relaxation(solve):
+    p = properties()
+    initial = np.full(5, .35)
+    thresholds = np.array([.06, .072, .14, .25, .30])
+    result = solve(np.zeros(60), np.full(60, 1e5), p,
+        pd.date_range("2020-01-01", periods=60), transpiration_data=np.full(60, 1e5),
+        initial_soil_moisture=initial, state_timing="end")
+    expected = thresholds + (initial-thresholds)*np.exp(-np.arange(1, 61)[:, None]/3.)
+    np.testing.assert_allclose(result, expected, rtol=0, atol=3e-16)
+    np.testing.assert_allclose(result.attrs["extraction_threshold"], thresholds, rtol=0, atol=1e-16)
+    np.testing.assert_array_equal(result.attrs["hydraulic_wilting_point"], p["wilting_point"])
+    assert np.all(result.iloc[-1, :3] < p["wilting_point"][:3])
+    assert np.all(result.values >= thresholds-1e-16)
+    b = result.attrs["water_balance"]
+    assert not np.any(b["storage_adjustment"]) and not np.any(b["actual_evaporation"])
+    changes = np.diff(np.vstack([initial, result.values]), axis=0)*p["layer_thickness"]
+    uptake = np.column_stack([b[f"transpiration_layer_{i}"] for i in range(1, 6)])
+    np.testing.assert_allclose(changes, -uptake, atol=3e-14, rtol=0)
+    np.testing.assert_allclose(changes.sum(axis=1), -b["actual_et"], atol=5e-14, rtol=0)
+
+
+def test_dry_initial_state_has_no_extraction_or_restoration(solve):
+    initial = np.array([.03, .036, .07, .125, .15])
+    result = solve([0.], [1000.], properties(), ["2020-01-01"],
+        transpiration_data=[1000.], initial_soil_moisture=initial, state_timing="end")
+    np.testing.assert_allclose(result.iloc[0], initial, rtol=0, atol=3e-17)
+    b = result.attrs["water_balance"]
+    assert b["actual_et"][0] == b["storage_adjustment"][0] == 0.
+    for layer in range(1, 6):
+        assert b[f"below_reference_days_layer_{layer}"][0] == 1.
+        assert b[f"below_wp_days_layer_{layer}"][0] == 1.
+
+
+def test_thresholds_do_not_block_redistribution_below_them(solve):
+    p = properties()
+    for key in ["layer_depth", "layer_thickness", "porosity", "wilting_point",
+                "available_water_capacity", "b_coefficient", "conductivity_sat"]:
+        p[key] = p[key][:2]
+    p.update(conductivity_sat=np.full(2, 100.), drainage_slope=0.)
+    initial = np.array([.05, .06])  # Both states are below their extraction thresholds.
+    result = solve([0.], [1000.], p, ["2020-01-01"], diff_factor=1000.,
+        transpiration_data=[1000.], initial_soil_moisture=initial, state_timing="end")
+    b = result.attrs["water_balance"]
+    assert b["interface_flux_1"][0] < 0.
+    assert result.iloc[0, 1] < initial[1]
+    assert b["actual_transpiration"][0] == 0.
+    assert abs(np.dot(result.iloc[0]-initial, p["layer_thickness"])) < 1e-14
+
+
+def test_unadjusted_reference_retains_wp_extraction_thresholds(solve):
+    p = properties()
+    p.pop("dry_reference_factor")
+    result = solve([0.], [1000.], p, ["2020-01-01"],
+        transpiration_data=[1000.], initial_soil_moisture=p["wilting_point"], state_timing="end")
+    np.testing.assert_array_equal(result.attrs["extraction_threshold"], p["wilting_point"])
+    assert result.attrs["water_balance"]["actual_transpiration"][0] == 0.
+
+
+def test_depth_rule_is_preserved_when_layers_cross_reference_bands(solve):
+    p = properties()
+    p.update(layer_depth=np.array([100., 200., 400., 1000.]),
+        layer_thickness=np.array([100., 100., 200., 600.]),
+        wilting_point=np.full(4, .2))
+    for key in ["porosity", "available_water_capacity", "b_coefficient", "conductivity_sat"]:
+        p[key] = p[key][:4]
+    result = solve([0.], [0.], p, ["2020-01-01"])
+    np.testing.assert_allclose(result.attrs["extraction_threshold"], [.08, .11, .17, .2], atol=1e-16)
+
+
+def test_zero_dry_reference_retains_the_third_layer_midpoint(solve):
+    p = properties()
+    p["dry_reference_factor"] = 0.
+    result = solve([0.], [0.], p, ["2020-01-01"])
+    np.testing.assert_array_equal(result.attrs["extraction_threshold"], [0., 0., .1, .25, .3])
+
+
+def test_repeated_and_disabled_initialization_preserve_revised_physics(solve):
+    """Combine the two branches against explicitly expanded daily forcing."""
+    p = properties()
+    p['conductivity_sat'] = np.full(5, 10.)
+    dates = pd.date_range('2014-01-01', '2016-12-31')
+    t = np.arange(len(dates))
+    rain = np.where(t % 11 == 0, 20., 0.)
+    et = 3. + np.sin(t / 50.)
+    tr = .8 * et
+    settings = dict(enabled=True, source_years=[2014, 2015], repeats=10)
+    compact = solve(rain, et, p, dates, transpiration_data=tr, spinup=settings)
+    expanded = [np.r_[np.tile(v[:730], 9), v] for v in (rain, et, tr)]
+    explicit = solve(expanded[0], expanded[1], p,
+        pd.date_range(end=dates[-1], periods=len(expanded[0])), transpiration_data=expanded[2])
+    np.testing.assert_allclose(compact, explicit.iloc[6570:], atol=1e-12, rtol=0)
+    np.testing.assert_allclose(compact.attrs['extraction_threshold'], [.06, .072, .14, .25, .30], atol=1e-16)
+    ordinary = solve(rain, et, p, dates, transpiration_data=tr)
+    disabled = solve(rain, et, p, dates, transpiration_data=tr, spinup={**settings, 'enabled': False})
+    np.testing.assert_array_equal(ordinary, disabled)
+    for key, values in compact.attrs['water_balance'].items():
+        np.testing.assert_allclose(values, explicit.attrs['water_balance'][key][6570:], atol=1e-10, rtol=0)
+        np.testing.assert_array_equal(ordinary.attrs['water_balance'][key], disabled.attrs['water_balance'][key])
diff --git a/tests/swb/test_multiyear_spinup.py b/tests/swb/test_multiyear_spinup.py
new file mode 100644
index 0000000..a4034fa
--- /dev/null
+++ b/tests/swb/test_multiyear_spinup.py
@@ -0,0 +1,209 @@
+"""
+Script: test_multiyear_spinup.py
+Objective: Verify compact repeated-year initialization against explicitly expanded histories.
+Author: Yi Yu (with assistance from Codex)
+Created: 2026-10-20
+Last updated: 2026-10-21
+Inputs: Deterministic daily forcing and physical five-layer columns.
+Outputs: Regression assertions for states, budgets, calendars and rejected configurations.
+Usage: pytest tests/swb/test_multiyear_spinup.py
+Dependencies: numpy, pandas, pytest, pysweb; numba for accelerated cases.
+"""
+import numpy as np
+import pandas as pd
+import pytest
+
+from pysweb.swb.spinup import source_block
+from pysweb.swb.solver import soil_water_balance_1d
+
+
+def properties():
+    bottoms = np.array([50., 150., 300., 600., 1000.])
+    return dict(layer_depth=bottoms, layer_thickness=np.diff(np.r_[0., bottoms]),
+                porosity=np.full(5, .45), wilting_point=np.full(5, .15),
+                available_water_capacity=np.full(5, .15), b_coefficient=np.full(5, 4.),
+                conductivity_sat=np.full(5, 10.), root_beta=.96)
+
+
+@pytest.mark.parametrize("backend", ["reference", "numba"])
+@pytest.mark.parametrize("years,repeats,timing", [
+    ([2014, 2015], 1, "start"), ([2014, 2015], 5, "start"),
+    ([2015, 2016], 3, "end"),
+])
+def test_compact_equals_explicit_history(backend, years, repeats, timing):
+    solve = soil_water_balance_1d
+    if backend == "numba":
+        pytest.importorskip("numba")
+        from pysweb.swb.accelerated import soil_water_balance_1d as solve
+    dates = pd.date_range(f"{years[0]}-01-01", f"{years[-1]+1}-02-17")
+    config = dict(enabled=True, source_years=years, repeats=repeats)
+    count, metadata = source_block(dates, config)
+    t = np.arange(len(dates))
+    rain = np.where(t % 9 == 0, 19., .1)
+    et = 2. + np.sin(t / 50.)
+    tr = .7 * et
+    compact = solve(rain, et, properties(), dates, transpiration_data=tr,
+                    spinup=config, state_timing=timing)
+    expanded = [np.r_[np.tile(v[:count], repeats-1), v] for v in (rain, et, tr)]
+    explicit = solve(expanded[0], expanded[1], properties(),
+                     pd.date_range(end=dates[-1], periods=len(expanded[0])),
+                     transpiration_data=expanded[2], state_timing=timing)
+    offset = count * (repeats-1)
+    np.testing.assert_allclose(compact, explicit.iloc[offset:], atol=1e-12, rtol=0)
+    np.testing.assert_allclose(compact.attrs["final_state"], explicit.attrs["final_state"], atol=1e-12)
+    for key, values in compact.attrs["water_balance"].items():
+        np.testing.assert_allclose(values, explicit.attrs["water_balance"][key][offset:], atol=1e-10)
+    assert compact.index.equals(dates)
+    assert metadata["initialization_days"] == count * repeats
+    assert compact.attrs["spinup"] == metadata
+
+
+@pytest.mark.parametrize("config", [
+    {"source_years": [2014, 2016], "repeats": 5},
+    {"source_years": [2014, 2015], "repeats": 0},
+    {"source_years": [2014, 2015], "repeats": True},
+    {"source_years": [2015, 2016], "repeats": 5},
+])
+def test_reject_incomplete_or_invalid_block(config):
+    with pytest.raises(ValueError):
+        source_block(pd.date_range("2014-01-01", "2016-12-31"), config)
+
+
+def test_missing_spinup_forcing_is_not_silently_skipped():
+    dates = pd.date_range("2014-01-01", "2016-01-02")
+    rain, et = np.ones(len(dates)), np.ones(len(dates))
+    rain[42] = np.nan
+    with pytest.raises(ValueError, match="invalid"):
+        soil_water_balance_1d(rain, et, properties(), dates,
+                             spinup=dict(enabled=True, source_years=[2014, 2015], repeats=5))
+
+
+def test_native_plan_excludes_source_years_from_scoring():
+    from pysweb.workflow_native import native_plan
+    config = dict(run=dict(start="2014-01-01", end="2025-12-31", state_timing="start",
+                           root="/tmp", name="fixture", holdout_start="2021-01-01"),
+                  calibration=dict(enabled=True, start="2016-01-01", end="2020-12-31"),
+                  spinup=dict(enabled=True, source_years=[2014, 2015], repeats=5,
+                              tolerance=.005, minimum_fraction=.99))
+    plan = native_plan(config, ["calibrate"])
+    assert plan["initialization"]["artificial_warmup_days"] == 4 * 730
+    assert plan["initialization"]["actual_warmup_days"] == 730
+    config["calibration"]["start"] = "2015-01-01"
+    with pytest.raises(ValueError, match="Scoring"):
+        native_plan(config, ["calibrate"])
+
+
+def test_failed_convergence_requires_explicit_record_action(tmp_path, monkeypatch):
+    import xarray as xr
+    import pysweb.workflow_protocol as protocol
+    import pysweb.swb.spinup as spinup_module
+    dates = pd.date_range("2014-01-01", "2015-12-31")
+    forcing = xr.Dataset({name: (("time", "lat", "lon"), np.ones((len(dates), 1, 1)))
+                          for name in ["effective_precip", "et", "E_t", "ndvi"]},
+                         coords={"time": dates, "lat": [0.], "lon": [0.]})
+    path = tmp_path / "forcing.nc"
+    forcing.to_netcdf(path)
+    soils = {"porosity": xr.DataArray(np.full((5, 1, 1), .45),
+             dims=("layer", "lat", "lon"), coords={"layer": np.arange(1, 6), "lat": [0.], "lon": [0.]})}
+    monkeypatch.setattr(protocol, "_load_soil_arrays", lambda *_: soils)
+    monkeypatch.setattr(protocol, "infer_layer_bottoms", lambda *_: [50, 150, 300, 600, 1000])
+    mask = soils["porosity"].isel(layer=0, drop=True).astype(bool)
+    monkeypatch.setattr(spinup_module, "verify_convergence", lambda *a, **k:
+                        (dict(passed=False, mode="wet_dry"), mask.astype(float)))
+    monkeypatch.setattr(protocol, "atomic_netcdf", lambda *a: None)
+    records = []
+    monkeypatch.setattr(protocol, "atomic_json", lambda p, r: records.append(dict(r)))
+    record = dict(scoring_start="2016-01-01", development_end="2020-12-31",
+                  actual_start="2014-01-01", simulation_start="2014-01-01",
+                  files={name: [str(path), name] for name in forcing})
+    config = dict(tolerance=.005, minimum_fraction=.99)
+    with pytest.raises(ValueError, match="convergence failed"):
+        protocol.convergence_gate(tmp_path, record, mask, {}, config, "default")
+    with pytest.warns(RuntimeWarning, match="did not converge"):
+        result = protocol.convergence_gate(tmp_path, record, mask, {},
+                    {**config, "convergence_action": "record"}, "default")
+    assert not result["passed"]
+    assert not records[-1]["passed"]
+
+
+@pytest.mark.parametrize("enabled", [True, False])
+def test_toml_to_native_plan_and_api_arguments(tmp_path, enabled):
+    from pysweb.workflow import load_config
+    from pysweb.workflow_native import native_plan
+    from pysweb.workflow_protocol import initialized_forcing
+    from pysweb.swb.run import _build_run_args
+    path = tmp_path / "experiment.toml"
+    calibration = '\n[calibration]\nenabled = true\nstart = "2016-01-01"\nend = "2020-12-31"\n' if enabled else ""
+    path.write_text('''[run]
+name = "fixture"
+root = "."
+extent = [146.2, -35.1, 146.3, -35.0]
+start = "2014-01-01"
+end = "2025-12-31"
+holdout_start = "2021-01-01"
+[spinup]
+enabled = true
+source_years = [2014, 2015]
+repeats = 10
+tolerance = 0.005
+minimum_fraction = 0.99
+convergence_action = "record"
+''' + calibration)
+    config = load_config(path)
+    plan = native_plan(config, ["run"])
+    assert plan["initialization"]["analysis_start"] == "2016-01-01"
+    assert plan["initialization"]["actual_warmup_days"] == 730
+    record = dict(simulation_start="2014-01-01", end="2025-12-31",
+                  initialization_mode="repeated_source_block", spinup=config["spinup"],
+                  files={key: [f"{key}.nc", key] for key in ["effective_precip", "et", "E_t", "ndvi"]})
+    arguments = initialized_forcing(record, tmp_path, workers=1, state_timing="start")
+    arguments.update(precip="precip.nc", output_dir=str(tmp_path))
+    args = _build_run_args(arguments)
+    assert args.spinup_source_years == [2014, 2015] and args.spinup_repeats == 10
+
+
+@pytest.mark.parametrize("backend", ["reference", "numba"])
+@pytest.mark.parametrize("configuration", [None, {}, {"enabled": False},
+    {"source_years": [2014, 2015], "repeats": 10},
+    {"enabled": False, "source_years": [2014, 2015], "repeats": 10}])
+def test_disabled_and_default_spinup_preserve_exact_unrepeated_run(backend, configuration):
+    solve = soil_water_balance_1d
+    if backend == "numba":
+        pytest.importorskip("numba")
+        from pysweb.swb.accelerated import soil_water_balance_1d as solve
+    dates = pd.date_range("2016-04-01", periods=9)
+    rain, et = np.arange(9.), np.full(9, 2.)
+    baseline = solve(rain, et, properties(), dates)
+    disabled = solve(rain, et, properties(), dates, spinup=configuration)
+    np.testing.assert_array_equal(disabled, baseline)
+    np.testing.assert_array_equal(disabled.attrs["initial_state"], baseline.attrs["initial_state"])
+    np.testing.assert_array_equal(disabled.attrs["final_state"], baseline.attrs["final_state"])
+    for key, values in baseline.attrs["water_balance"].items():
+        np.testing.assert_array_equal(disabled.attrs["water_balance"][key], values)
+    assert "spinup" not in disabled.attrs
+
+
+@pytest.mark.parametrize("flag", ["false", "true", 0, 1, None])
+def test_spinup_switch_rejects_non_boolean_values(flag):
+    from pysweb.swb.spinup import enabled_spinup
+    with pytest.raises(ValueError, match="boolean"):
+        enabled_spinup({"enabled": flag})
+
+
+@pytest.mark.parametrize("table", ["", "[spinup]\n", "[spinup]\nenabled = false\n",
+    "[spinup]\nsource_years = [2014, 2015]\nrepeats = 10\n"])
+def test_disabled_toml_needs_no_source_years_or_development_boundary(tmp_path, table):
+    from pysweb.workflow import load_config
+    from pysweb.workflow_native import native_plan
+    path = tmp_path / "disabled.toml"
+    path.write_text('''[run]
+name = "fixture"
+root = "."
+extent = [146.2, -35.1, 146.3, -35.0]
+start = "2016-04-01"
+end = "2016-04-09"
+''' + table)
+    plan = native_plan(load_config(path), ["run"])
+    assert plan["initialization"]["repeated_spinup_enabled"] is False
+    assert plan["initialization"]["artificial_warmup_days"] == 0
+    assert plan["initialization"]["analysis_start"] == "2016-04-01"
diff --git a/tests/swb/test_references.py b/tests/swb/test_references.py
index 7c7a840..14d5b55 100644
--- a/tests/swb/test_references.py
+++ b/tests/swb/test_references.py
@@ -3,7 +3,7 @@
 Objective: Verify temporal-before-spatial GSSM quantiles, support integrity and artifact identity.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-10-18
-Last updated: 2026-10-18
+Last updated: 2026-10-21
 Inputs: Synthetic geographic cells with anticorrelated temporal moisture extremes.
 Outputs: Exact reference order, coverage and tamper-rejection assertions.
 Usage: python -m pytest tests/swb/test_references.py
@@ -61,12 +61,13 @@ def test_missing_soil_or_incomplete_reference_cannot_shrink_frozen_support():
         derive_reference(gssm, mask, hyd, site="fixture", sources=["synthetic"], hydraulic_provenance="fixture")
 
 
-def test_daily_model_reuses_unchanged_observation_reference(tmp_path):
+@pytest.mark.parametrize("previous", ["swb-conservative-2p-v1", "swb-daily-conservative-2p-v2"])
+def test_daily_model_reuses_unchanged_observation_reference(tmp_path, previous):
     from pysweb.swb.formulation import canonical_id
     gssm, mask, hyd = inputs()
     record = derive_reference(gssm, mask, hyd, site="fixture", sources=["synthetic"],
                               hydraulic_provenance="fixture")
-    record["formulation_version"] = "swb-conservative-2p-v1"
+    record["formulation_version"] = previous
     record["artifact_id"] = canonical_id(record)
     path = tmp_path / "reference.json"
     path.write_text(json.dumps(record))
diff --git a/tests/test_protocols.py b/tests/test_protocols.py
index b3f0b6f..9b737ac 100644
--- a/tests/test_protocols.py
+++ b/tests/test_protocols.py
@@ -3,7 +3,7 @@
 Objective: Verify holdout isolation, donor restrictions, spin-up and evaluation contracts.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-10-17
-Last updated: 2026-10-18
+Last updated: 2026-10-21
 Inputs: Small deterministic synthetic scenes, forcing and observation arrays.
 Outputs: Scientific protocol regression assertions.
 Usage: python -m pytest tests/test_protocols.py
@@ -380,7 +380,8 @@ def test_screened_scene_workflow_offline(tmp_path):
             reconstruction=dict(acquisition_max_gap_days=64,boundary_carry=True),workers=1)
 
 
-def test_configured_spinup_products_and_endpoint_operator(tmp_path):
+@pytest.mark.parametrize("spinup_enabled", [True, False])
+def test_configured_spinup_products_and_endpoint_operator(tmp_path, spinup_enabled):
     """Run real orchestration, convergence, fitting, reproduction and extraction."""
     from tests.io.test_native_grids import raster
     from rasterio import Affine
@@ -427,18 +428,23 @@ def test_configured_spinup_products_and_endpoint_operator(tmp_path):
             no_polish=True,serial_polish=True,solver_backend=SOLVER_BACKEND,seed=0,diff_bounds=[1000.,1000.],
             beta_bounds=[.96,.96]),
         model=dict(solver_backend=SOLVER_BACKEND),forcing=dict(peer_runs=[]),
-        spinup=dict(source_year=2014,repeats=2,tolerance=.005,minimum_fraction=.99,solver_backend=SOLVER_BACKEND))
+        spinup=dict(enabled=spinup_enabled,source_year=2014,repeats=2,tolerance=.005,
+                    minimum_fraction=.99,solver_backend=SOLVER_BACKEND))
     run_native(config,['benchmark'])
     assert not (base/'calibration.csv').exists()
     benchmark=json.loads((base/'calibration-benchmark.json').read_text())
-    assert benchmark['cells']==4 and benchmark['simulated_days']==1095
+    assert benchmark['cells']==4 and benchmark['simulated_days']==(1095 if spinup_enabled else 365)
     run_native(config,['calibrate','run'])
     acceptance=json.loads((base/'calibration-acceptance.json').read_text())
     assert acceptance['passed'] and acceptance['same_pair_count']
     for folder in (coarse,fine):
         for label in ('default','calibrated'):
-            gate=json.loads((folder/'initialized'/f'{label}_convergence.json').read_text())
-            assert gate['passed'] and gate['initialization_days']==730
+            path=folder/'initialized'/f'{label}_convergence.json'
+            if spinup_enabled:
+                gate=json.loads(path.read_text())
+                assert gate['passed'] and gate['initialization_days']==730
+            else:
+                assert not path.exists()
     product=base/'products/SWEB_RZSM.nc'
     longitude,latitude=Transformer.from_crs(projected.rio.crs,4326,always_xy=True).transform(
         projected.x.values[0],projected.y.values[0])
@@ -446,7 +452,11 @@ def test_configured_spinup_products_and_endpoint_operator(tmp_path):
     assert len(table)==len(dates) and table.continuous_trajectory_valid.all()
     with xr.open_dataset(product) as ds:
         assert ds.sizes['time']==len(dates)
-        assert ds.initial_history_valid_default.astype(bool).all()
+        if spinup_enabled:
+            assert ds.initial_history_valid_default.astype(bool).all()
+        else:
+            assert 'initial_history_valid_default' not in ds
+            np.testing.assert_allclose(ds.rzsm_default.isel(time=0), .32, atol=1e-7)
         last=ds.rzsm_calibrated.isel(time=-1,y=0,x=0).values
         final=np.array([ds[f'final_sm_layer_{i}_calibrated'].isel(y=0,x=0).item() for i in (1,2)])
         assert table.storage_dailymean_mm.iloc[-1]==pytest.approx(((last+final)/2)@np.array([20.,150.]))