From b25be94c858090d92ee1e5593a43b813bb120655 Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Fri, 28 Aug 2026 21:37:57 +1000 Subject: [PATCH] Add an optional coarse-root pool to FF16 and TF24 Both models carried one root pool, fine roots proportional to leaf area, so below-ground structure did not scale with tree size. Coarse roots are now the below-ground continuation of the sapwood cylinder, m_cr = a_cr1 * m_s, which makes the structural root:shoot ratio size-invariant and keeps the growth allocation closed-form. a_cr1 defaults to 0, so no model's output and no scientific version moves. expand_state() gains a mass_coarse_root column. Closes #349 Co-authored-by: Claude Opus 5 (1M context) --- NEWS.md | 14 ++ R/RcppR6.R | 2 +- inst/RcppR6_classes.yml | 6 + inst/docs/FF16/FF16-eqns.csv | 3 +- inst/docs/FF16/FF16-params-core.csv | 5 +- inst/include/plant/RcppR6_post.hpp | 18 ++ .../plant/models/ff16_production_kernel.h | 27 ++- inst/include/plant/models/ff16_strategy.h | 38 +++- inst/include/plant/models/tf24_strategy.h | 40 +++- src/ff16_strategy.cpp | 72 +++++-- src/strategy_expand.cpp | 11 +- src/tf24_strategy.cpp | 71 +++++-- tests/testthat/_snaps/model-version.md | 12 +- tests/testthat/test-coarse-roots.R | 181 ++++++++++++++++++ tests/testthat/test-expand-state.R | 24 ++- tests/testthat/test-ff16-ad-kernel.R | 12 +- .../testthat/test-ff16-resident-coupling-ad.R | 6 +- .../test-strategy-ff16-reference-comparison.R | 7 +- tests/testthat/test-strategy-ff16.R | 3 + tests/testthat/test-strategy-tf24.R | 3 + tests/testthat/test-strategy-tf24f.R | 3 + 21 files changed, 487 insertions(+), 71 deletions(-) create mode 100644 tests/testthat/test-coarse-roots.R diff --git a/NEWS.md b/NEWS.md index 92188367f..ffdbca0d1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -440,6 +440,20 @@ were not previously recorded here: ### New features +* **FF16 and TF24 gained a coarse-root pool, off by default (#349).** Both models carried one root pool, `mass_root = a_r1 * A_l` — fine roots, with no height scaling. Coarse (structural) roots are now modelled as the below-ground continuation of the sapwood cylinder: + + ``` + m_cr = a_cr1 * m_s = a_cr1 * rho * theta * eta_c * A_l * H + ``` + + so they scale with height and the structural root:shoot ratio is size-invariant — the property fine roots cannot supply. They respire at `r_cr` and turn over at `k_cr`, both defaulting to the sapwood rates because coarse roots are woody, and they enter `mass_live` / `mass_total` but not `mass_above_ground`. `*_expand_state()` gains a `mass_coarse_root` column. + + Because `d m_cr / d A_l = a_cr1 * d m_s / d A_l`, the growth allocation stays closed-form: the new term reuses the sapwood term the bark share already shares, and `darea_leaf_dmass_live` gains one addend rather than needing a rederivation. + + **`a_cr1` defaults to 0, so nothing moves and no scientific version does either.** Every new term is appended rather than interleaved, so at the default the sums are bit-identical (`x + 0.0 == x`, and the association order of the existing terms is untouched) — `FF16@v1`, `TF24@v9` and `TF24f@v9.1` all stand, and FF16 still agrees with the published 2012 reference model. The model-version drift guard fires on the three added parameters; its snapshot was re-blessed after confirming the diff is three added names and no changed default. + + ⚠️ **Coarse roots are deliberately outside TF24's water uptake.** The depth-distributed uptake reads `a_r1` directly: absorbing surface is a fine-root property, and `test-coarse-roots.R` asserts every water-side quantity is bit-identical with the pool on. + * **The NSC storage pool is bounded by the shape of its own flow (`TF24@v9`, `TF24f@v9.1`).** `dS/dt` was `net_flux > 0 ? net_flux : floor_gate * net_flux` -- one signed flux with a gate applied in one direction. It is now a charge and diff --git a/R/RcppR6.R b/R/RcppR6.R index 95df43e9e..e8a731471 100644 --- a/R/RcppR6.R +++ b/R/RcppR6.R @@ -1,6 +1,6 @@ ## Generated by RcppR6: do not edit by hand ## Version: 0.2.4 -## Hash: 1a68109cddbb80490e6a4d4d186276fd +## Hash: e8df3bbeda019508b5c7ac6024ad00cc ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index a848d841d..37f2c879d 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -1150,10 +1150,12 @@ FF16_Pars: - a_l1: double - a_l2: double - a_r1: double + - a_cr1: double - a_b1: double - r_s: double - r_b: double - r_r: double + - r_cr: double - r_l: double - a_y: double - a_bio: double @@ -1161,6 +1163,7 @@ FF16_Pars: - k_b: double - k_s: double - k_r: double + - k_cr: double - a_p1: double - a_p2: double - a_f3: double @@ -1319,10 +1322,12 @@ TF24_Pars: - a_l1: double - a_l2: double - a_r1: double + - a_cr1: double - a_b1: double - r_s: double - r_b: double - r_r: double + - r_cr: double - r_l: double - a_y: double - a_bio: double @@ -1330,6 +1335,7 @@ TF24_Pars: - k_b: double - k_s: double - k_r: double + - k_cr: double - a_p1: double - a_p2: double - a_f3: double diff --git a/inst/docs/FF16/FF16-eqns.csv b/inst/docs/FF16/FF16-eqns.csv index 0cee2377d..5c231f4b9 100644 --- a/inst/docs/FF16/FF16-eqns.csv +++ b/inst/docs/FF16/FF16-eqns.csv @@ -7,5 +7,6 @@ "Leaf mass","$M_{\rm l} = \phi \, A_{\rm l}$","$\frac{{\rm d}M_{\rm l}}{{\rm d}A_{\rm l}} = \phi$","$\frac{{\rm d}M_{\rm l}}{{\rm d}t} = \frac{{\rm d}M_{\rm l}}{{\rm d}A_{\rm l}} \, \frac{{\rm d}A_{\rm l}}{{\rm d}t}$" "Sapwood mass","$M_{\rm s} = \rho \, \theta \, \eta_{\rm c} \, A_{\rm l} \, H$","$\frac{{\rm d}M_{\rm s}}{{\rm d}A_{\rm l}} = \rho \, \theta \, \eta_{\rm c} \, \big(H + A_{\rm l} \, \frac{{\rm d}H}{{\rm d}A_{\rm l}} \big)$","$\frac{{\rm d}M_{\rm s}}{{\rm d}t} = \frac{{\rm d}M_{\rm s}}{{\rm d}A_{\rm l}} \, \frac{{\rm d}A_{\rm l}}{{\rm d}t}$" "Bark mass","$M_{\rm b} = \alpha_{{\rm b1}} \, \rho \, \theta \, \eta_{\rm c} \, A_{\rm l} \, H$","$\frac{{\rm d}M_{\rm b}}{{\rm d}A_{\rm l}} = \alpha_{{\rm b1}} \, \rho \, \theta \, \eta_{\rm c}\big(H + A_{\rm l} \, \frac{{\rm d}H}{{\rm d}A_{\rm l}} \big)$","$\frac{{\rm d}M_{\rm b}}{{\rm d}t} = \frac{{\rm d}M_{\rm b}}{{\rm d}A_{\rm l}} \, \frac{{\rm d}A_{\rm l}}{{\rm d}t}$" -"Root mass","$M_{\rm r} = \alpha_{{\rm r1}} \, A_{\rm l}$","$\frac{{\rm d}M_{\rm r}}{{\rm d}A_{\rm l}} = \alpha_{{\rm r1}}$","$\frac{{\rm d}M_{\rm r}}{{\rm d}t} = \frac{{\rm d}M_{\rm r}}{{\rm d}A_{\rm l}} \, \frac{{\rm d}A_{\rm l}}{{\rm d}t}$" +"Fine-root mass","$M_{\rm r} = \alpha_{{\rm r1}} \, A_{\rm l}$","$\frac{{\rm d}M_{\rm r}}{{\rm d}A_{\rm l}} = \alpha_{{\rm r1}}$","$\frac{{\rm d}M_{\rm r}}{{\rm d}t} = \frac{{\rm d}M_{\rm r}}{{\rm d}A_{\rm l}} \, \frac{{\rm d}A_{\rm l}}{{\rm d}t}$" +"Coarse-root mass","$M_{\rm cr} = \alpha_{{\rm cr1}} \, M_{\rm s}$","$\frac{{\rm d}M_{\rm cr}}{{\rm d}A_{\rm l}} = \alpha_{{\rm cr1}} \, \frac{{\rm d}M_{\rm s}}{{\rm d}A_{\rm l}}$","$\frac{{\rm d}M_{\rm cr}}{{\rm d}t} = \frac{{\rm d}M_{\rm cr}}{{\rm d}A_{\rm l}} \, \frac{{\rm d}A_{\rm l}}{{\rm d}t}$" diff --git a/inst/docs/FF16/FF16-params-core.csv b/inst/docs/FF16/FF16-params-core.csv index 0e1e02510..1a659ddce 100644 --- a/inst/docs/FF16/FF16-params-core.csv +++ b/inst/docs/FF16/FF16-params-core.csv @@ -6,7 +6,8 @@ Wood density,$\rho$,kg m$^{-3}$,rho Sapwood area per unit leaf area,$\theta$,,theta Height of plant with leaf area of 1m$^2$,$\alpha_{\rm l1}$,m,a_l1 Exponent of relationship between height and leaf area,$\alpha_{\rm l2}$,,a_l2 -Root mass per unit leaf area,$\alpha_{\rm r1}$,kg m$^{-2}$,a_r1 +Fine-root mass per unit leaf area,$\alpha_{\rm r1}$,kg m$^{-2}$,a_r1 +Coarse-root mass per unit sapwood mass,$\alpha_{\rm cr1}$,,a_cr1 Ratio of bark area to sapwood area,$\alpha_{\rm b1}$,,a_b1 **Production**,,, Leaf photosynthesis per area,$\alpha_{\rm p1}$,mol yr$^{-1}$ m$^{-2}$,a_p1 @@ -17,10 +18,12 @@ Leaf respiration per mass,$r_{\rm l}$,mol yr$^{-1}$ kg$^{-1}$,r_l Fine-root respiration per mass,$r_{\rm r}$,mol yr$^{-1}$ kg$^{-1}$,r_r Sapwood respiration per mass,$r_{\rm s}$,mol yr$^{-1}$ kg$^{-1}$,r_s Bark respiration per mass,$r_{\rm b}$,mol yr$^{-1}$ kg$^{-1}$,r_b +Coarse-root respiration per mass,$r_{\rm cr}$,mol yr$^{-1}$ kg$^{-1}$,r_cr Turnover rate for leaves,$k_{\rm l}$,yr$^{-1}$,k_l Turnover rate for sapwood,$k_{\rm s}$,yr$^{-1}$,k_s Turnover rate for bark,$k_{\rm b}$,yr$^{-1}$,k_b Turnover rate for fine roots,$k_{\rm r}$,yr$^{-1}$,k_r +Turnover rate for coarse roots,$k_{\rm cr}$,yr$^{-1}$,k_cr **Fecundity**,,, Seed mass,$\omega$,kg,omega Height at maturation,$H_{{\rm mat}}$,m,hmat diff --git a/inst/include/plant/RcppR6_post.hpp b/inst/include/plant/RcppR6_post.hpp index f53cc4d70..374fb8da3 100644 --- a/inst/include/plant/RcppR6_post.hpp +++ b/inst/include/plant/RcppR6_post.hpp @@ -1105,10 +1105,12 @@ template <> inline SEXP wrap(const plant::FF16_Pars& x) { ret["a_l1"] = Rcpp::wrap(x.a_l1); ret["a_l2"] = Rcpp::wrap(x.a_l2); ret["a_r1"] = Rcpp::wrap(x.a_r1); + ret["a_cr1"] = Rcpp::wrap(x.a_cr1); ret["a_b1"] = Rcpp::wrap(x.a_b1); ret["r_s"] = Rcpp::wrap(x.r_s); ret["r_b"] = Rcpp::wrap(x.r_b); ret["r_r"] = Rcpp::wrap(x.r_r); + ret["r_cr"] = Rcpp::wrap(x.r_cr); ret["r_l"] = Rcpp::wrap(x.r_l); ret["a_y"] = Rcpp::wrap(x.a_y); ret["a_bio"] = Rcpp::wrap(x.a_bio); @@ -1116,6 +1118,7 @@ template <> inline SEXP wrap(const plant::FF16_Pars& x) { ret["k_b"] = Rcpp::wrap(x.k_b); ret["k_s"] = Rcpp::wrap(x.k_s); ret["k_r"] = Rcpp::wrap(x.k_r); + ret["k_cr"] = Rcpp::wrap(x.k_cr); ret["a_p1"] = Rcpp::wrap(x.a_p1); ret["a_p2"] = Rcpp::wrap(x.a_p2); ret["a_f3"] = Rcpp::wrap(x.a_f3); @@ -1158,6 +1161,8 @@ template <> inline plant::FF16_Pars as(SEXP x) { ret.a_l2 = Rcpp::as(xl["a_l2"]); // ret.a_r1 = Rcpp::as(xl["a_r1"]); ret.a_r1 = Rcpp::as(xl["a_r1"]); + // ret.a_cr1 = Rcpp::as(xl["a_cr1"]); + ret.a_cr1 = Rcpp::as(xl["a_cr1"]); // ret.a_b1 = Rcpp::as(xl["a_b1"]); ret.a_b1 = Rcpp::as(xl["a_b1"]); // ret.r_s = Rcpp::as(xl["r_s"]); @@ -1166,6 +1171,8 @@ template <> inline plant::FF16_Pars as(SEXP x) { ret.r_b = Rcpp::as(xl["r_b"]); // ret.r_r = Rcpp::as(xl["r_r"]); ret.r_r = Rcpp::as(xl["r_r"]); + // ret.r_cr = Rcpp::as(xl["r_cr"]); + ret.r_cr = Rcpp::as(xl["r_cr"]); // ret.r_l = Rcpp::as(xl["r_l"]); ret.r_l = Rcpp::as(xl["r_l"]); // ret.a_y = Rcpp::as(xl["a_y"]); @@ -1180,6 +1187,8 @@ template <> inline plant::FF16_Pars as(SEXP x) { ret.k_s = Rcpp::as(xl["k_s"]); // ret.k_r = Rcpp::as(xl["k_r"]); ret.k_r = Rcpp::as(xl["k_r"]); + // ret.k_cr = Rcpp::as(xl["k_cr"]); + ret.k_cr = Rcpp::as(xl["k_cr"]); // ret.a_p1 = Rcpp::as(xl["a_p1"]); ret.a_p1 = Rcpp::as(xl["a_p1"]); // ret.a_p2 = Rcpp::as(xl["a_p2"]); @@ -1349,10 +1358,12 @@ template <> inline SEXP wrap(const plant::TF24_Pars& x) { ret["a_l1"] = Rcpp::wrap(x.a_l1); ret["a_l2"] = Rcpp::wrap(x.a_l2); ret["a_r1"] = Rcpp::wrap(x.a_r1); + ret["a_cr1"] = Rcpp::wrap(x.a_cr1); ret["a_b1"] = Rcpp::wrap(x.a_b1); ret["r_s"] = Rcpp::wrap(x.r_s); ret["r_b"] = Rcpp::wrap(x.r_b); ret["r_r"] = Rcpp::wrap(x.r_r); + ret["r_cr"] = Rcpp::wrap(x.r_cr); ret["r_l"] = Rcpp::wrap(x.r_l); ret["a_y"] = Rcpp::wrap(x.a_y); ret["a_bio"] = Rcpp::wrap(x.a_bio); @@ -1360,6 +1371,7 @@ template <> inline SEXP wrap(const plant::TF24_Pars& x) { ret["k_b"] = Rcpp::wrap(x.k_b); ret["k_s"] = Rcpp::wrap(x.k_s); ret["k_r"] = Rcpp::wrap(x.k_r); + ret["k_cr"] = Rcpp::wrap(x.k_cr); ret["a_p1"] = Rcpp::wrap(x.a_p1); ret["a_p2"] = Rcpp::wrap(x.a_p2); ret["a_f3"] = Rcpp::wrap(x.a_f3); @@ -1432,6 +1444,8 @@ template <> inline plant::TF24_Pars as(SEXP x) { ret.a_l2 = Rcpp::as(xl["a_l2"]); // ret.a_r1 = Rcpp::as(xl["a_r1"]); ret.a_r1 = Rcpp::as(xl["a_r1"]); + // ret.a_cr1 = Rcpp::as(xl["a_cr1"]); + ret.a_cr1 = Rcpp::as(xl["a_cr1"]); // ret.a_b1 = Rcpp::as(xl["a_b1"]); ret.a_b1 = Rcpp::as(xl["a_b1"]); // ret.r_s = Rcpp::as(xl["r_s"]); @@ -1440,6 +1454,8 @@ template <> inline plant::TF24_Pars as(SEXP x) { ret.r_b = Rcpp::as(xl["r_b"]); // ret.r_r = Rcpp::as(xl["r_r"]); ret.r_r = Rcpp::as(xl["r_r"]); + // ret.r_cr = Rcpp::as(xl["r_cr"]); + ret.r_cr = Rcpp::as(xl["r_cr"]); // ret.r_l = Rcpp::as(xl["r_l"]); ret.r_l = Rcpp::as(xl["r_l"]); // ret.a_y = Rcpp::as(xl["a_y"]); @@ -1454,6 +1470,8 @@ template <> inline plant::TF24_Pars as(SEXP x) { ret.k_s = Rcpp::as(xl["k_s"]); // ret.k_r = Rcpp::as(xl["k_r"]); ret.k_r = Rcpp::as(xl["k_r"]); + // ret.k_cr = Rcpp::as(xl["k_cr"]); + ret.k_cr = Rcpp::as(xl["k_cr"]); // ret.a_p1 = Rcpp::as(xl["a_p1"]); ret.a_p1 = Rcpp::as(xl["a_p1"]); // ret.a_p2 = Rcpp::as(xl["a_p2"]); diff --git a/inst/include/plant/models/ff16_production_kernel.h b/inst/include/plant/models/ff16_production_kernel.h index 702956c33..fd0d06983 100644 --- a/inst/include/plant/models/ff16_production_kernel.h +++ b/inst/include/plant/models/ff16_production_kernel.h @@ -32,17 +32,23 @@ S ff16_assimilation_leaf(S a_p1, S a_p2, S x) { } // [eqn 13] Total maintenance respiration (linear in the mass cascade). +// The coarse-root term is appended last so that at the default a_cr1 = 0 (hence +// mass_coarse_root == 0) the sum is bit-identical to the pre-#349 expression. template S ff16_respiration(S mass_leaf, S mass_sapwood, S mass_bark, S mass_root, - S r_l, S r_s, S r_b, S r_r) { - return r_l * mass_leaf + r_b * mass_bark + r_s * mass_sapwood + r_r * mass_root; + S mass_coarse_root, + S r_l, S r_s, S r_b, S r_r, S r_cr) { + return r_l * mass_leaf + r_b * mass_bark + r_s * mass_sapwood + r_r * mass_root + + r_cr * mass_coarse_root; } // [eqn 14] Total turnover. template S ff16_turnover(S mass_leaf, S mass_bark, S mass_sapwood, S mass_root, - S k_l, S k_b, S k_s, S k_r) { - return k_l * mass_leaf + k_b * mass_bark + k_s * mass_sapwood + k_r * mass_root; + S mass_coarse_root, + S k_l, S k_b, S k_s, S k_r, S k_cr) { + return k_l * mass_leaf + k_b * mass_bark + k_s * mass_sapwood + k_r * mass_root + + k_cr * mass_coarse_root; } // [eqn 15] Net production from assimilation/respiration/turnover. @@ -57,10 +63,10 @@ S ff16_net_production_A(S a_bio, S a_y, S assimilation, S respiration, S turnove // functions directly with pars.* members. template struct FF16ProdPars { - S lma, rho, theta, a_b1, a_r1, eta_c; + S lma, rho, theta, a_b1, a_r1, a_cr1, eta_c; S a_p1, a_p2; - S r_l, r_s, r_b, r_r; - S k_l, k_b, k_s, k_r; + S r_l, r_s, r_b, r_r, r_cr; + S k_l, k_b, k_s, k_r, k_cr; S a_bio, a_y; }; @@ -82,11 +88,14 @@ S ff16_net_from_components(const FF16ProdPars& p, S height, S area_leaf, const S area_bark = p.a_b1 * area_leaf * p.theta; const S mass_bark = area_bark * height * p.eta_c * p.rho; const S mass_root = p.a_r1 * area_leaf; + const S mass_coarse_root = p.a_cr1 * mass_sapwood; const S respiration = ff16_respiration(mass_leaf, mass_sapwood, mass_bark, mass_root, - p.r_l, p.r_s, p.r_b, p.r_r); + mass_coarse_root, + p.r_l, p.r_s, p.r_b, p.r_r, p.r_cr); const S turnover = ff16_turnover(mass_leaf, mass_bark, mass_sapwood, mass_root, - p.k_l, p.k_b, p.k_s, p.k_r); + mass_coarse_root, + p.k_l, p.k_b, p.k_s, p.k_r, p.k_cr); return ff16_net_production_A(p.a_bio, p.a_y, assimilation, respiration, turnover); } diff --git a/inst/include/plant/models/ff16_strategy.h b/inst/include/plant/models/ff16_strategy.h index b7e988565..2394d82ff 100644 --- a/inst/include/plant/models/ff16_strategy.h +++ b/inst/include/plant/models/ff16_strategy.h @@ -32,6 +32,12 @@ struct FF16_Pars { double a_l2 = 0.306; // dimensionless scaling of height with leaf area // Root mass per leaf area double a_r1 = 0.07; //[kg / m] + // Coarse (structural) root mass per unit sapwood mass. Coarse roots are + // treated as the below-ground continuation of the sapwood cylinder, so they + // inherit its height scaling and the structural root:shoot ratio is + // size-invariant. Zero by default: the pool is absent unless switched on, so + // FF16 reproduces the published model (Falster et al. 2016). + double a_cr1 = 0.0; // [dimensionless] // Ratio of bark area : sapwood area double a_b1 = 0.17; // [dimensionless] @@ -54,6 +60,10 @@ struct FF16_Pars { // -- r_b is defined (new) as 2*r_s, whereas the paper assumes a // fixed multiplication by 2) double r_b = 2.0 * r_s; + // Coarse-root respiration per mass [mol CO2 / yr / kg] + // Coarse roots are woody, so this defaults to the sapwood rate. As with r_b, + // this is a default, not a link: setting r_s afterwards does not move it. + double r_cr = r_s; // Carbon conversion parameter double a_y = 0.7; // Constant converting assimilated CO2 to dry mass [kg / mol] @@ -67,6 +77,9 @@ struct FF16_Pars { double k_s = 0.2; // Root turnover [/yr] double k_r = 1.0; + // Coarse-root turnover [/yr] + // Woody, so this defaults to the sapwood rate rather than the fine-root one. + double k_cr = k_s; // Parameters of the hyperbola for annual LRC double a_p1 = 151.177775377968; // [mol CO2 / yr / m2] double a_p2 = 0.204716166503633; // [dimensionless] @@ -192,15 +205,21 @@ class FF16_Strategy: public Strategy { // [eqn 7] Mass of (fine) roots double mass_root(double area_leaf) const; + // [eqn 7b] Mass of coarse (structural) roots, as a fraction of sapwood mass + // (see FF16_Pars::a_cr1). Zero unless a_cr1 is set. + double mass_coarse_root(double mass_sapwood) const; + // [eqn 8] Total Mass double mass_live(double mass_leaf, double mass_bark, - double mass_sapwood, double mass_root) const; + double mass_sapwood, double mass_root, + double mass_coarse_root) const; double mass_total(double mass_leaf, double mass_bark, double mass_sapwood, - double mass_heartwood, double mass_root) const; + double mass_heartwood, double mass_root, + double mass_coarse_root) const; // Above-ground mass = leaf + all stem components (bark + sapwood + - // heartwood); excludes roots. + // heartwood); excludes both fine and coarse roots. double mass_above_ground(double mass_leaf, double mass_bark, double mass_sapwood, double mass_heartwood) const; @@ -257,20 +276,24 @@ class FF16_Strategy: public Strategy { // [eqn 13] Total maintenance respiration double respiration(double mass_leaf, double mass_sapwood, - double mass_bark, double mass_root) const; + double mass_bark, double mass_root, + double mass_coarse_root) const; double respiration_leaf(double mass) const; double respiration_bark(double mass) const; double respiration_sapwood(double mass) const; double respiration_root(double mass) const; + double respiration_coarse_root(double mass) const; // [eqn 14] Total turnover double turnover(double mass_leaf, double mass_bark, - double mass_sapwood, double mass_root) const; + double mass_sapwood, double mass_root, + double mass_coarse_root) const; double turnover_leaf(double mass) const; double turnover_bark(double mass) const; double turnover_sapwood(double mass) const; double turnover_root(double mass) const; + double turnover_coarse_root(double mass) const; // [eqn 15] Net production double net_mass_production_dt_A(double assimilation, double respiration, @@ -326,6 +349,9 @@ class FF16_Strategy: public Strategy { double dmass_bark_darea_leaf(double area_leaf) const; // Mass of root needed for new unit area leaf, d m_r / d a_l double dmass_root_darea_leaf(double area_leaf) const; + // Mass of coarse root needed for new unit area leaf, d m_cr / d a_l. + // == a_cr1 * dmass_sapwood_darea_leaf, so it carries the same height scaling. + double dmass_coarse_root_darea_leaf(double area_leaf) const; // Growth rate of basal diameter_stem per unit stem area double ddiameter_stem_darea_stem(double area_stem) const; // Growth rate of components per unit time: @@ -337,6 +363,8 @@ class FF16_Strategy: public Strategy { double diameter_stem_dt(double area_stem, double area_stem_dt) const; double mass_root_dt(double area_leaf, double area_leaf_dt) const; + double mass_coarse_root_dt(double area_leaf, + double area_leaf_dt) const; double mass_live_dt(double fraction_allocation_reproduction, double net_mass_production_dt) const; double mass_total_dt(double fraction_allocation_reproduction, diff --git a/inst/include/plant/models/tf24_strategy.h b/inst/include/plant/models/tf24_strategy.h index e9baa984d..ba44b973d 100644 --- a/inst/include/plant/models/tf24_strategy.h +++ b/inst/include/plant/models/tf24_strategy.h @@ -29,19 +29,30 @@ struct TF24_Pars { double theta = 1.0/4669; // Sapwood area per leaf area [dimensionless] double a_l1 = 5.44; // height with 1m2 leaf [m] double a_l2 = 0.306; // scaling of height with leaf area - double a_r1 = 0.07; // Root mass per leaf area [kg / m] + double a_r1 = 0.07; // Fine-root mass per leaf area [kg / m] + // Coarse (structural) root mass per unit sapwood mass. Coarse roots are + // treated as the below-ground continuation of the sapwood cylinder, so they + // inherit its height scaling and the structural root:shoot ratio is + // size-invariant. Zero by default, so the pool is absent unless switched on. + // Only a_r1 feeds the depth-distributed water uptake -- that is a fine-root + // process, and coarse roots must stay out of it (see #349). + double a_cr1 = 0.0; // [dimensionless] double a_b1 = 0.17; // Ratio of bark area : sapwood area // * Production double r_s = 4012.0 / 608.0; // Sapwood respiration per stem mass double r_b = 2.0 * r_s; // Bark respiration (assumed 2 x sapwood) - double r_r = 217.0; // Root respiration per mass + double r_r = 217.0; // Fine-root respiration per mass + // Coarse roots are woody, so this defaults to the sapwood rate. As with r_b, + // a default and not a link: setting r_s afterwards does not move it. + double r_cr = r_s; // Coarse-root respiration per mass double r_l = 39.27 / 0.1978791; // Leaf dark respiration per leaf mass double a_y = 0.7; // Carbon conversion parameter double a_bio = 2.45e-2; // CO2 -> dry mass [kg / mol] double k_l = 0.4565855; // Leaf turnover [/yr] double k_b = 0.2; // Bark turnover [/yr] double k_s = 0.2; // Sapwood turnover [/yr] - double k_r = 1.0; // Root turnover [/yr] + double k_r = 1.0; // Fine-root turnover [/yr] + double k_cr = k_s; // Coarse-root turnover [/yr] (woody) double a_p1 = 151.177775377968; // LRC hyperbola [mol CO2 / yr / m2] double a_p2 = 0.204716166503633; // LRC hyperbola shape // * Seed production @@ -378,12 +389,18 @@ class TF24_Strategy: public Strategy { // [eqn 7] Mass of (fine) roots double mass_root(double area_leaf) const; + // [eqn 7b] Mass of coarse (structural) roots, as a fraction of sapwood mass + // (see TF24_Pars::a_cr1). Zero unless a_cr1 is set. + double mass_coarse_root(double mass_sapwood) const; + // [eqn 8] Total Mass double mass_live(double mass_leaf, double mass_bark, - double mass_sapwood, double mass_root) const; + double mass_sapwood, double mass_root, + double mass_coarse_root) const; double mass_total(double mass_leaf, double mass_bark, double mass_sapwood, - double mass_heartwood, double mass_root) const; + double mass_heartwood, double mass_root, + double mass_coarse_root) const; // Above-ground mass = leaf + all stem components (bark + sapwood + // heartwood); excludes roots. @@ -407,20 +424,24 @@ class TF24_Strategy: public Strategy { // [eqn 13] Total maintenance respiration double respiration(double mass_leaf, double mass_sapwood, - double mass_bark, double mass_root) const; + double mass_bark, double mass_root, + double mass_coarse_root) const; double respiration_leaf(double mass) const; double respiration_bark(double mass) const; double respiration_sapwood(double mass) const; double respiration_root(double mass) const; + double respiration_coarse_root(double mass) const; // [eqn 14] Total turnover double turnover(double mass_leaf, double mass_bark, - double mass_sapwood, double mass_root) const; + double mass_sapwood, double mass_root, + double mass_coarse_root) const; double turnover_leaf(double mass) const; double turnover_bark(double mass) const; double turnover_sapwood(double mass) const; double turnover_root(double mass) const; + double turnover_coarse_root(double mass) const; // [eqn 15] Net production double net_mass_production_dt_A(double assimilation, double respiration, @@ -467,6 +488,9 @@ class TF24_Strategy: public Strategy { double dmass_bark_darea_leaf(double area_leaf) const; // Mass of root needed for new unit area leaf, d m_r / d a_l double dmass_root_darea_leaf(double area_leaf) const; + // Mass of coarse root needed for new unit area leaf, d m_cr / d a_l. + // == a_cr1 * dmass_sapwood_darea_leaf, so it carries the same height scaling. + double dmass_coarse_root_darea_leaf(double area_leaf) const; // Growth rate of basal diameter_stem per unit stem area double ddiameter_stem_darea_stem(double area_stem) const; // Growth rate of components per unit time: @@ -478,6 +502,8 @@ class TF24_Strategy: public Strategy { double diameter_stem_dt(double area_stem, double area_stem_dt) const; double mass_root_dt(double area_leaf, double area_leaf_dt) const; + double mass_coarse_root_dt(double area_leaf, + double area_leaf_dt) const; double mass_live_dt(double fraction_allocation_reproduction, double net_mass_production_dt) const; double mass_total_dt(double fraction_allocation_reproduction, diff --git a/src/ff16_strategy.cpp b/src/ff16_strategy.cpp index 40fafe3d0..98faa87a3 100644 --- a/src/ff16_strategy.cpp +++ b/src/ff16_strategy.cpp @@ -63,16 +63,26 @@ double FF16_Strategy::mass_root(double area_leaf) const { return pars.a_r1 * area_leaf; } +// [eqn 7b] Mass of coarse (structural) roots +double FF16_Strategy::mass_coarse_root(double mass_sapwood) const { + return pars.a_cr1 * mass_sapwood; +} + // [eqn 8] Total mass +// The coarse-root term is appended, not interleaved, so that with the default +// a_cr1 = 0 the sum is bit-identical to the pre-#349 expression (x + 0.0 == x +// exactly, and the association order of the other terms is untouched). double FF16_Strategy::mass_live(double mass_leaf, double mass_bark, - double mass_sapwood, double mass_root) const { - return mass_leaf + mass_sapwood + mass_bark + mass_root; + double mass_sapwood, double mass_root, + double mass_coarse_root) const { + return mass_leaf + mass_sapwood + mass_bark + mass_root + mass_coarse_root; } double FF16_Strategy::mass_total(double mass_leaf, double mass_bark, double mass_sapwood, double mass_heartwood, - double mass_root) const { - return mass_leaf + mass_bark + mass_sapwood + mass_heartwood + mass_root; + double mass_root, double mass_coarse_root) const { + return mass_leaf + mass_bark + mass_sapwood + mass_heartwood + mass_root + + mass_coarse_root; } double FF16_Strategy::mass_above_ground(double mass_leaf, double mass_bark, @@ -210,10 +220,12 @@ double FF16_Strategy::assimilation_leaf(double x) const { // [eqn 13] Total maintenance respiration // NOTE: In contrast with Falster ref model, we do not normalise by pars.a_y*pars.a_bio. double FF16_Strategy::respiration(double mass_leaf, double mass_sapwood, - double mass_bark, double mass_root) const { + double mass_bark, double mass_root, + double mass_coarse_root) const { // Single source: scalar-templated kernel (#472 scope B, Milestone A). return ff16_respiration(mass_leaf, mass_sapwood, mass_bark, mass_root, - pars.r_l, pars.r_s, pars.r_b, pars.r_r); + mass_coarse_root, + pars.r_l, pars.r_s, pars.r_b, pars.r_r, pars.r_cr); } double FF16_Strategy::respiration_leaf(double mass) const { @@ -232,12 +244,18 @@ double FF16_Strategy::respiration_root(double mass) const { return pars.r_r * mass; } +double FF16_Strategy::respiration_coarse_root(double mass) const { + return pars.r_cr * mass; +} + // [eqn 14] Total turnover double FF16_Strategy::turnover(double mass_leaf, double mass_bark, - double mass_sapwood, double mass_root) const { + double mass_sapwood, double mass_root, + double mass_coarse_root) const { // Single source: scalar-templated kernel (#472 scope B, Milestone A). return ff16_turnover(mass_leaf, mass_bark, mass_sapwood, mass_root, - pars.k_l, pars.k_b, pars.k_s, pars.k_r); + mass_coarse_root, + pars.k_l, pars.k_b, pars.k_s, pars.k_r, pars.k_cr); } double FF16_Strategy::turnover_leaf(double mass) const { @@ -256,6 +274,10 @@ double FF16_Strategy::turnover_root(double mass) const { return pars.k_r * mass; } +double FF16_Strategy::turnover_coarse_root(double mass) const { + return pars.k_cr * mass; +} + // [eqn 15] Net production // // NOTE: Translation of variable names from the Falster 2011. Everything @@ -292,12 +314,15 @@ double FF16_Strategy::net_mass_production_dt(const FF16_Environment& environment const double area_bark_ = area_bark(area_leaf_); const double mass_bark_ = mass_bark(area_bark_, height); const double mass_root_ = mass_root(area_leaf_); + const double mass_coarse_root_ = mass_coarse_root(mass_sapwood_); const double assimilation_ = assimilation(environment, height, area_leaf_, height_inverse); const double respiration_ = - respiration(mass_leaf_, mass_sapwood_, mass_bark_, mass_root_); + respiration(mass_leaf_, mass_sapwood_, mass_bark_, mass_root_, + mass_coarse_root_); const double turnover_ = - turnover(mass_leaf_, mass_bark_, mass_sapwood_, mass_root_); + turnover(mass_leaf_, mass_bark_, mass_sapwood_, mass_root_, + mass_coarse_root_); return net_mass_production_dt_A(assimilation_, respiration_, turnover_); } @@ -328,10 +353,14 @@ double FF16_Strategy::darea_leaf_dmass_live(double area_leaf, // so compute the shared pow(area_leaf, pars.a_l2) term once rather than twice. const double dmass_sapwood_darea_leaf_ = dmass_sapwood_darea_leaf(area_leaf, area_leaf_pow_a_l2); + // dmass_coarse_root_darea_leaf == a_cr1 * dmass_sapwood_darea_leaf, so it + // reuses the same shared term as bark. Appended last so that at the default + // a_cr1 = 0 the denominator is bit-identical to the pre-#349 sum. return 1.0/( dmass_leaf_darea_leaf(area_leaf) + dmass_sapwood_darea_leaf_ + pars.a_b1 * dmass_sapwood_darea_leaf_ - + dmass_root_darea_leaf(area_leaf)); + + dmass_root_darea_leaf(area_leaf) + + pars.a_cr1 * dmass_sapwood_darea_leaf_); } double FF16_Strategy::dheight_darea_leaf(double area_leaf) const { @@ -371,6 +400,11 @@ double FF16_Strategy::dmass_root_darea_leaf(double /* area_leaf */) const { return pars.a_r1; } +// Mass of coarse root needed for new unit area leaf, d m_cr / d a_l +double FF16_Strategy::dmass_coarse_root_darea_leaf(double area_leaf) const { + return pars.a_cr1 * dmass_sapwood_darea_leaf(area_leaf); +} + // Growth rate of basal diameter_stem per unit time double FF16_Strategy::ddiameter_stem_darea_stem(double area_stem) const { return pow(M_PI * area_stem, -0.5); @@ -411,6 +445,12 @@ double FF16_Strategy::mass_root_dt(double area_leaf, return area_leaf_dt * dmass_root_darea_leaf(area_leaf); } +// Growth rate of coarse-root mass per unit time +double FF16_Strategy::mass_coarse_root_dt(double area_leaf, + double area_leaf_dt) const { + return area_leaf_dt * dmass_coarse_root_darea_leaf(area_leaf); +} + double FF16_Strategy::mass_live_dt(double fraction_allocation_reproduction, double net_mass_production_dt) const { return (1 - fraction_allocation_reproduction) * net_mass_production_dt; @@ -431,8 +471,10 @@ double FF16_Strategy::mass_above_ground_dt(double area_leaf, double area_leaf_dt) const { const double mass_root_dt = area_leaf_dt * dmass_root_darea_leaf(area_leaf); + const double mass_coarse_root_dt = + area_leaf_dt * dmass_coarse_root_darea_leaf(area_leaf); return mass_total_dt(fraction_allocation_reproduction, net_mass_production_dt, - mass_heartwood_dt) - mass_root_dt; + mass_heartwood_dt) - mass_root_dt - mass_coarse_root_dt; } double FF16_Strategy::mass_heartwood_dt(double mass_sapwood) const { @@ -442,10 +484,12 @@ double FF16_Strategy::mass_heartwood_dt(double mass_sapwood) const { double FF16_Strategy::mass_live_given_height(double height) const { double area_leaf_ = area_leaf(height); + const double mass_sapwood_ = mass_sapwood(area_sapwood(area_leaf_), height); return mass_leaf(area_leaf_) + mass_bark(area_bark(area_leaf_), height) + - mass_sapwood(area_sapwood(area_leaf_), height) + - mass_root(area_leaf_); + mass_sapwood_ + + mass_root(area_leaf_) + + mass_coarse_root(mass_sapwood_); } double FF16_Strategy::height_given_mass_leaf(double mass_leaf) const { diff --git a/src/strategy_expand.cpp b/src/strategy_expand.cpp index 348130315..5e4bad49d 100644 --- a/src/strategy_expand.cpp +++ b/src/strategy_expand.cpp @@ -24,8 +24,8 @@ Rcpp::List strategy_expand_allometry_impl(S s, } Rcpp::NumericVector area_leaf(n), mass_leaf(n), area_sapwood(n), mass_sapwood(n), area_bark(n), mass_bark(n), area_stem(n), - diameter_stem(n), mass_root(n), mass_live(n), mass_total(n), - mass_above_ground(n); + diameter_stem(n), mass_root(n), mass_coarse_root(n), mass_live(n), + mass_total(n), mass_above_ground(n); for (R_xlen_t i = 0; i < n; ++i) { const double h = height[i]; @@ -37,6 +37,7 @@ Rcpp::List strategy_expand_allometry_impl(S s, const double mb = s.mass_bark(ab, h); const double ast = s.area_stem(ab, as, area_heartwood[i]); const double mr = s.mass_root(al); + const double mcr = s.mass_coarse_root(ms); const double mh = mass_heartwood[i]; area_leaf[i] = al; @@ -48,8 +49,9 @@ Rcpp::List strategy_expand_allometry_impl(S s, area_stem[i] = ast; diameter_stem[i] = s.diameter_stem(ast); mass_root[i] = mr; - mass_live[i] = s.mass_live(ml, mb, ms, mr); - mass_total[i] = s.mass_total(ml, mb, ms, mh, mr); + mass_coarse_root[i] = mcr; + mass_live[i] = s.mass_live(ml, mb, ms, mr, mcr); + mass_total[i] = s.mass_total(ml, mb, ms, mh, mr, mcr); mass_above_ground[i] = s.mass_above_ground(ml, mb, ms, mh); } @@ -63,6 +65,7 @@ Rcpp::List strategy_expand_allometry_impl(S s, Rcpp::_["area_stem"] = area_stem, Rcpp::_["diameter_stem"] = diameter_stem, Rcpp::_["mass_root"] = mass_root, + Rcpp::_["mass_coarse_root"] = mass_coarse_root, Rcpp::_["mass_live"] = mass_live, Rcpp::_["mass_total"] = mass_total, Rcpp::_["mass_above_ground"] = mass_above_ground); diff --git a/src/tf24_strategy.cpp b/src/tf24_strategy.cpp index 97684ba24..dbb842809 100644 --- a/src/tf24_strategy.cpp +++ b/src/tf24_strategy.cpp @@ -123,16 +123,26 @@ double TF24_Strategy::mass_root(double area_leaf) const { return pars.a_r1 * area_leaf; } +// [eqn 7b] Mass of coarse (structural) roots +double TF24_Strategy::mass_coarse_root(double mass_sapwood) const { + return pars.a_cr1 * mass_sapwood; +} + // [eqn 8] Total mass +// The coarse-root term is appended, not interleaved, so that with the default +// a_cr1 = 0 the sum is bit-identical to the pre-#349 expression (x + 0.0 == x +// exactly, and the association order of the other terms is untouched). double TF24_Strategy::mass_live(double mass_leaf, double mass_bark, - double mass_sapwood, double mass_root) const { - return mass_leaf + mass_sapwood + mass_bark + mass_root; + double mass_sapwood, double mass_root, + double mass_coarse_root) const { + return mass_leaf + mass_sapwood + mass_bark + mass_root + mass_coarse_root; } double TF24_Strategy::mass_total(double mass_leaf, double mass_bark, double mass_sapwood, double mass_heartwood, - double mass_root) const { - return mass_leaf + mass_bark + mass_sapwood + mass_heartwood + mass_root; + double mass_root, double mass_coarse_root) const { + return mass_leaf + mass_bark + mass_sapwood + mass_heartwood + mass_root + + mass_coarse_root; } double TF24_Strategy::mass_above_ground(double mass_leaf, double mass_bark, @@ -326,11 +336,13 @@ double TF24_Strategy::assimilation_leaf(double x) const { // [eqn 13] Total maintenance respiration // NOTE: In contrast with Falster ref model, we do not normalise by pars.a_y*pars.a_bio. double TF24_Strategy::respiration(double mass_leaf, double mass_sapwood, - double mass_bark, double mass_root) const { + double mass_bark, double mass_root, + double mass_coarse_root) const { return respiration_leaf(mass_leaf) + respiration_bark(mass_bark) + respiration_sapwood(mass_sapwood) + - respiration_root(mass_root); + respiration_root(mass_root) + + respiration_coarse_root(mass_coarse_root); } double TF24_Strategy::respiration_leaf(double mass) const { @@ -349,13 +361,19 @@ double TF24_Strategy::respiration_root(double mass) const { return pars.r_r * mass; } +double TF24_Strategy::respiration_coarse_root(double mass) const { + return pars.r_cr * mass; +} + // [eqn 14] Total turnover double TF24_Strategy::turnover(double mass_leaf, double mass_bark, - double mass_sapwood, double mass_root) const { + double mass_sapwood, double mass_root, + double mass_coarse_root) const { return turnover_leaf(mass_leaf) + turnover_bark(mass_bark) + turnover_sapwood(mass_sapwood) + - turnover_root(mass_root); + turnover_root(mass_root) + + turnover_coarse_root(mass_coarse_root); } double TF24_Strategy::turnover_leaf(double mass) const { @@ -374,6 +392,10 @@ double TF24_Strategy::turnover_root(double mass) const { return pars.k_r * mass; } +double TF24_Strategy::turnover_coarse_root(double mass) const { + return pars.k_cr * mass; +} + // [eqn 15] Net production // // NOTE: Translation of variable names from the Falster 2011. Everything @@ -609,10 +631,13 @@ double TF24_Strategy::net_mass_production_dt(const TF24_Environment& environment // converts to canopy area, then years, then mols const double assimilation_ = carbon_profit_ * area_leaf_* 60*60*12*365/1e6; // const double assimilation_ = assimilation(environment, height, area_leaf_); + const double mass_coarse_root_ = mass_coarse_root(mass_sapwood_); const double respiration_ = - respiration(mass_leaf_, mass_sapwood_, mass_bark_, mass_root_); + respiration(mass_leaf_, mass_sapwood_, mass_bark_, mass_root_, + mass_coarse_root_); const double turnover_ = - turnover(mass_leaf_, mass_bark_, mass_sapwood_, mass_root_); + turnover(mass_leaf_, mass_bark_, mass_sapwood_, mass_root_, + mass_coarse_root_); return net_mass_production_dt_A(assimilation_, respiration_, turnover_); } @@ -662,10 +687,13 @@ double TF24_Strategy::fecundity_dt(double net_mass_production_dt, } double TF24_Strategy::darea_leaf_dmass_live(double area_leaf) const { + // The coarse-root term is appended last so that at the default a_cr1 = 0 the + // denominator is bit-identical to the pre-#349 sum. return 1.0/( dmass_leaf_darea_leaf(area_leaf) + dmass_sapwood_darea_leaf(area_leaf) + dmass_bark_darea_leaf(area_leaf) - + dmass_root_darea_leaf(area_leaf)); + + dmass_root_darea_leaf(area_leaf) + + dmass_coarse_root_darea_leaf(area_leaf)); } double TF24_Strategy::dheight_darea_leaf(double area_leaf) const { @@ -692,6 +720,11 @@ double TF24_Strategy::dmass_root_darea_leaf(double /* area_leaf */) const { return pars.a_r1; } +// Mass of coarse root needed for new unit area leaf, d m_cr / d a_l +double TF24_Strategy::dmass_coarse_root_darea_leaf(double area_leaf) const { + return pars.a_cr1 * dmass_sapwood_darea_leaf(area_leaf); +} + // Growth rate of basal diameter_stem per unit time double TF24_Strategy::ddiameter_stem_darea_stem(double area_stem) const { return pow(M_PI * area_stem, -0.5); @@ -732,6 +765,12 @@ double TF24_Strategy::mass_root_dt(double area_leaf, return area_leaf_dt * dmass_root_darea_leaf(area_leaf); } +// Growth rate of coarse-root mass per unit time +double TF24_Strategy::mass_coarse_root_dt(double area_leaf, + double area_leaf_dt) const { + return area_leaf_dt * dmass_coarse_root_darea_leaf(area_leaf); +} + double TF24_Strategy::mass_live_dt(double fraction_allocation_reproduction, double net_mass_production_dt) const { return (1 - fraction_allocation_reproduction) * net_mass_production_dt; @@ -752,8 +791,10 @@ double TF24_Strategy::mass_above_ground_dt(double area_leaf, double area_leaf_dt) const { const double mass_root_dt = area_leaf_dt * dmass_root_darea_leaf(area_leaf); + const double mass_coarse_root_dt = + area_leaf_dt * dmass_coarse_root_darea_leaf(area_leaf); return mass_total_dt(fraction_allocation_reproduction, net_mass_production_dt, - mass_heartwood_dt) - mass_root_dt; + mass_heartwood_dt) - mass_root_dt - mass_coarse_root_dt; } double TF24_Strategy::mass_heartwood_dt(double mass_sapwood) const { @@ -763,10 +804,12 @@ double TF24_Strategy::mass_heartwood_dt(double mass_sapwood) const { double TF24_Strategy::mass_live_given_height(double height) const { double area_leaf_ = area_leaf(height); + const double mass_sapwood_ = mass_sapwood(area_sapwood(area_leaf_), height); return mass_leaf(area_leaf_) + mass_bark(area_bark(area_leaf_), height) + - mass_sapwood(area_sapwood(area_leaf_), height) + - mass_root(area_leaf_); + mass_sapwood_ + + mass_root(area_leaf_) + + mass_coarse_root(mass_sapwood_); } double TF24_Strategy::height_given_mass_leaf(double mass_leaf) const { diff --git a/tests/testthat/_snaps/model-version.md b/tests/testthat/_snaps/model-version.md index b184c1816..84c256a88 100644 --- a/tests/testthat/_snaps/model-version.md +++ b/tests/testthat/_snaps/model-version.md @@ -6,10 +6,10 @@ "names": { "type": "character", "attributes": {}, - "value": ["birth_rate_x", "birth_rate_y", "control.GSS_tol_abs", "control.ci_abs_tol", "control.ci_niter", "control.fixed_time_step", "control.function_integration_rule", "control.node_density_in_birth_date", "control.node_gradient_direction", "control.node_gradient_eps", "control.node_gradient_richardson", "control.node_gradient_richardson_depth", "control.ode_a_dydt", "control.ode_a_y", "control.ode_step_size_initial", "control.ode_step_size_max", "control.ode_step_size_min", "control.ode_tol_abs", "control.ode_tol_rel", "control.offspring_production_iterations", "control.offspring_production_tol", "control.ppa_layer_optical_depth", "control.ppa_layer_smoothing", "control.save_RK45_cache", "control.schedule_eps", "control.schedule_nsteps", "control.schedule_verbose", "control.shading_model", "control.vulnerability_curve_ncontrol", "is_variable_birth_rate", "model_id", "pars.S_D", "pars.a_b1", "pars.a_bio", "pars.a_d0", "pars.a_dG1", "pars.a_dG2", "pars.a_f1", "pars.a_f2", "pars.a_f3", "pars.a_l1", "pars.a_l2", "pars.a_p1", "pars.a_p2", "pars.a_r1", "pars.a_y", "pars.d_I", "pars.eta", "pars.hmat", "pars.k_I", "pars.k_b", "pars.k_l", "pars.k_r", "pars.k_s", "pars.lma", "pars.omega", "pars.r_b", "pars.r_l", "pars.r_r", "pars.r_s", "pars.recruitment_decay", "pars.rho", "pars.theta"] + "value": ["birth_rate_x", "birth_rate_y", "control.GSS_tol_abs", "control.ci_abs_tol", "control.ci_niter", "control.fixed_time_step", "control.function_integration_rule", "control.node_density_in_birth_date", "control.node_gradient_direction", "control.node_gradient_eps", "control.node_gradient_richardson", "control.node_gradient_richardson_depth", "control.ode_a_dydt", "control.ode_a_y", "control.ode_step_size_initial", "control.ode_step_size_max", "control.ode_step_size_min", "control.ode_tol_abs", "control.ode_tol_rel", "control.offspring_production_iterations", "control.offspring_production_tol", "control.ppa_layer_optical_depth", "control.ppa_layer_smoothing", "control.save_RK45_cache", "control.schedule_eps", "control.schedule_nsteps", "control.schedule_verbose", "control.shading_model", "control.vulnerability_curve_ncontrol", "is_variable_birth_rate", "model_id", "pars.S_D", "pars.a_b1", "pars.a_bio", "pars.a_cr1", "pars.a_d0", "pars.a_dG1", "pars.a_dG2", "pars.a_f1", "pars.a_f2", "pars.a_f3", "pars.a_l1", "pars.a_l2", "pars.a_p1", "pars.a_p2", "pars.a_r1", "pars.a_y", "pars.d_I", "pars.eta", "pars.hmat", "pars.k_I", "pars.k_b", "pars.k_cr", "pars.k_l", "pars.k_r", "pars.k_s", "pars.lma", "pars.omega", "pars.r_b", "pars.r_cr", "pars.r_l", "pars.r_r", "pars.r_s", "pars.recruitment_decay", "pars.rho", "pars.theta"] } }, - "value": ["", "1", "0.001", "0.001", "1000", "0", "21", "FALSE", "-1", "9.9999999999999995e-07", "FALSE", "4", "0", "1", "9.9999999999999995e-07", "5", "9.9999999999999995e-07", "0.0001", "0.0001", "1000", "1e-08", "0.5", "0.29999999999999999", "FALSE", "0.02", "20", "FALSE", "", "100", "FALSE", "FF16@v1", "0.25", "0.17000000000000001", "0.024500000000000001", "0.10000000000000001", "5.5", "20", "1", "50", "0.00011400000000000001", "5.4400000000000004", "0.30599999999999999", "151.177775377968", "0.204716166503633", "0.070000000000000007", "0.69999999999999996", "0.01", "12", "16.595869100000002", "0.5", "0.20000000000000001", "0.45658549999999998", "1", "0.20000000000000001", "0.1978791", "3.8000000000000002e-05", "13.197368421052632", "198.45451086041933", "217", "6.5986842105263159", "0", "608", "0.00021417862497322766"] + "value": ["", "1", "0.001", "0.001", "1000", "0", "21", "FALSE", "-1", "9.9999999999999995e-07", "FALSE", "4", "0", "1", "9.9999999999999995e-07", "5", "9.9999999999999995e-07", "0.0001", "0.0001", "1000", "1e-08", "0.5", "0.29999999999999999", "FALSE", "0.02", "20", "FALSE", "", "100", "FALSE", "FF16@v1", "0.25", "0.17000000000000001", "0.024500000000000001", "0", "0.10000000000000001", "5.5", "20", "1", "50", "0.00011400000000000001", "5.4400000000000004", "0.30599999999999999", "151.177775377968", "0.204716166503633", "0.070000000000000007", "0.69999999999999996", "0.01", "12", "16.595869100000002", "0.5", "0.20000000000000001", "0.20000000000000001", "0.45658549999999998", "1", "0.20000000000000001", "0.1978791", "3.8000000000000002e-05", "13.197368421052632", "6.5986842105263159", "198.45451086041933", "217", "6.5986842105263159", "0", "608", "0.00021417862497322766"] } # scientific surface of K93 matches its declared version @@ -34,10 +34,10 @@ "names": { "type": "character", "attributes": {}, - "value": ["birth_rate_x", "birth_rate_y", "control.GSS_tol_abs", "control.ci_abs_tol", "control.ci_niter", "control.fixed_time_step", "control.function_integration_rule", "control.node_density_in_birth_date", "control.node_gradient_direction", "control.node_gradient_eps", "control.node_gradient_richardson", "control.node_gradient_richardson_depth", "control.ode_a_dydt", "control.ode_a_y", "control.ode_step_size_initial", "control.ode_step_size_max", "control.ode_step_size_min", "control.ode_tol_abs", "control.ode_tol_rel", "control.offspring_production_iterations", "control.offspring_production_tol", "control.ppa_layer_optical_depth", "control.ppa_layer_smoothing", "control.save_RK45_cache", "control.schedule_eps", "control.schedule_nsteps", "control.schedule_verbose", "control.shading_model", "control.vulnerability_curve_ncontrol", "is_variable_birth_rate", "model_id", "pars.K_s", "pars.S_D", "pars.TF24_beta2", "pars.TF24_cost_scale", "pars.TF24_floor_lambda_o", "pars.a", "pars.a_b1", "pars.a_bio", "pars.a_d0", "pars.a_dG1", "pars.a_dG2", "pars.a_f1", "pars.a_f2", "pars.a_f3", "pars.a_l1", "pars.a_l2", "pars.a_p1", "pars.a_p2", "pars.a_r1", "pars.a_st1", "pars.a_st2", "pars.a_st3", "pars.a_y", "pars.beta1", "pars.curv_fact_colim", "pars.curv_fact_elec_trans", "pars.d", "pars.d_I", "pars.dmass_dN", "pars.eta", "pars.hmat", "pars.jmax_25", "pars.k_I", "pars.k_b", "pars.k_l", "pars.k_r", "pars.k_s", "pars.lma", "pars.nmass_b", "pars.nmass_l", "pars.nmass_r", "pars.nmass_s", "pars.omega", "pars.psi_crit", "pars.r_b", "pars.r_l", "pars.r_r", "pars.r_s", "pars.recruitment_decay", "pars.rho", "pars.root_b", "pars.root_c", "pars.root_depth_shape_eta", "pars.root_psi_crit", "pars.rooting_depth_max", "pars.stem_P50", "pars.stem_b", "pars.stem_c", "pars.theta", "pars.use_energy_balance", "pars.var_sapwood_volume_cost", "pars.vcmax_25"] + "value": ["birth_rate_x", "birth_rate_y", "control.GSS_tol_abs", "control.ci_abs_tol", "control.ci_niter", "control.fixed_time_step", "control.function_integration_rule", "control.node_density_in_birth_date", "control.node_gradient_direction", "control.node_gradient_eps", "control.node_gradient_richardson", "control.node_gradient_richardson_depth", "control.ode_a_dydt", "control.ode_a_y", "control.ode_step_size_initial", "control.ode_step_size_max", "control.ode_step_size_min", "control.ode_tol_abs", "control.ode_tol_rel", "control.offspring_production_iterations", "control.offspring_production_tol", "control.ppa_layer_optical_depth", "control.ppa_layer_smoothing", "control.save_RK45_cache", "control.schedule_eps", "control.schedule_nsteps", "control.schedule_verbose", "control.shading_model", "control.vulnerability_curve_ncontrol", "is_variable_birth_rate", "model_id", "pars.K_s", "pars.S_D", "pars.TF24_beta2", "pars.TF24_cost_scale", "pars.TF24_floor_lambda_o", "pars.a", "pars.a_b1", "pars.a_bio", "pars.a_cr1", "pars.a_d0", "pars.a_dG1", "pars.a_dG2", "pars.a_f1", "pars.a_f2", "pars.a_f3", "pars.a_l1", "pars.a_l2", "pars.a_p1", "pars.a_p2", "pars.a_r1", "pars.a_st1", "pars.a_st2", "pars.a_st3", "pars.a_y", "pars.beta1", "pars.curv_fact_colim", "pars.curv_fact_elec_trans", "pars.d", "pars.d_I", "pars.dmass_dN", "pars.eta", "pars.hmat", "pars.jmax_25", "pars.k_I", "pars.k_b", "pars.k_cr", "pars.k_l", "pars.k_r", "pars.k_s", "pars.lma", "pars.nmass_b", "pars.nmass_l", "pars.nmass_r", "pars.nmass_s", "pars.omega", "pars.psi_crit", "pars.r_b", "pars.r_cr", "pars.r_l", "pars.r_r", "pars.r_s", "pars.recruitment_decay", "pars.rho", "pars.root_b", "pars.root_c", "pars.root_depth_shape_eta", "pars.root_psi_crit", "pars.rooting_depth_max", "pars.stem_P50", "pars.stem_b", "pars.stem_c", "pars.theta", "pars.use_energy_balance", "pars.var_sapwood_volume_cost", "pars.vcmax_25"] } }, - "value": ["", "1", "0.001", "0.001", "1000", "0", "21", "FALSE", "-1", "9.9999999999999995e-07", "FALSE", "4", "0", "1", "9.9999999999999995e-07", "5", "9.9999999999999995e-07", "0.0001", "0.0001", "1000", "1e-08", "0.5", "0.29999999999999999", "FALSE", "0.02", "20", "FALSE", "", "100", "FALSE", "TF24@v9", "1", "0.25", "1.5", "7.5", "0", "0.29999999999999999", "0.17000000000000001", "0.024500000000000001", "0.10000000000000001", "5.5", "20", "1", "50", "0.00011400000000000001", "5.4400000000000004", "0.30599999999999999", "151.177775377968", "0.204716166503633", "0.070000000000000007", "0.10000000000000001", "0.10000000000000001", "0.80000000000000004", "0.69999999999999996", "20000", "0.98999999999999999", "0.69999999999999996", "0.050000000000000003", "0.01", "0", "12", "16.595869100000002", "157.44", "0.5", "0.20000000000000001", "0.45658549999999998", "1", "0.20000000000000001", "0.1978791", "0.0033999999999999998", "0.012999999999999999", "0.0033500000000000001", "0.00198", "3.8000000000000002e-05", "7.0854928681315172", "13.197368421052632", "198.45451086041933", "217", "6.5986842105263159", "0", "608", "3.8982450000000002", "2.6801469999999998", "0.20000000000000001", "5.8702825428827037", "1.5", "1.8500000000000001", "2.5894373183081751", "1.0899851679673995", "0.00021417862497322766", "0", "1", "96"] + "value": ["", "1", "0.001", "0.001", "1000", "0", "21", "FALSE", "-1", "9.9999999999999995e-07", "FALSE", "4", "0", "1", "9.9999999999999995e-07", "5", "9.9999999999999995e-07", "0.0001", "0.0001", "1000", "1e-08", "0.5", "0.29999999999999999", "FALSE", "0.02", "20", "FALSE", "", "100", "FALSE", "TF24@v9", "1", "0.25", "1.5", "7.5", "0", "0.29999999999999999", "0.17000000000000001", "0.024500000000000001", "0", "0.10000000000000001", "5.5", "20", "1", "50", "0.00011400000000000001", "5.4400000000000004", "0.30599999999999999", "151.177775377968", "0.204716166503633", "0.070000000000000007", "0.10000000000000001", "0.10000000000000001", "0.80000000000000004", "0.69999999999999996", "20000", "0.98999999999999999", "0.69999999999999996", "0.050000000000000003", "0.01", "0", "12", "16.595869100000002", "157.44", "0.5", "0.20000000000000001", "0.20000000000000001", "0.45658549999999998", "1", "0.20000000000000001", "0.1978791", "0.0033999999999999998", "0.012999999999999999", "0.0033500000000000001", "0.00198", "3.8000000000000002e-05", "7.0854928681315172", "13.197368421052632", "6.5986842105263159", "198.45451086041933", "217", "6.5986842105263159", "0", "608", "3.8982450000000002", "2.6801469999999998", "0.20000000000000001", "5.8702825428827037", "1.5", "1.8500000000000001", "2.5894373183081751", "1.0899851679673995", "0.00021417862497322766", "0", "1", "96"] } # scientific surface of TF24f matches its declared version @@ -48,9 +48,9 @@ "names": { "type": "character", "attributes": {}, - "value": ["birth_rate_x", "birth_rate_y", "control.GSS_tol_abs", "control.ci_abs_tol", "control.ci_niter", "control.fixed_time_step", "control.function_integration_rule", "control.node_density_in_birth_date", "control.node_gradient_direction", "control.node_gradient_eps", "control.node_gradient_richardson", "control.node_gradient_richardson_depth", "control.ode_a_dydt", "control.ode_a_y", "control.ode_step_size_initial", "control.ode_step_size_max", "control.ode_step_size_min", "control.ode_tol_abs", "control.ode_tol_rel", "control.offspring_production_iterations", "control.offspring_production_tol", "control.ppa_layer_optical_depth", "control.ppa_layer_smoothing", "control.save_RK45_cache", "control.schedule_eps", "control.schedule_nsteps", "control.schedule_verbose", "control.shading_model", "control.vulnerability_curve_ncontrol", "is_variable_birth_rate", "model_id", "pars.K_s", "pars.S_D", "pars.TF24_beta2", "pars.TF24_cost_scale", "pars.TF24_floor_lambda_o", "pars.a", "pars.a_b1", "pars.a_bio", "pars.a_d0", "pars.a_dG1", "pars.a_dG2", "pars.a_f1", "pars.a_f2", "pars.a_f3", "pars.a_l1", "pars.a_l2", "pars.a_p1", "pars.a_p2", "pars.a_r1", "pars.a_st1", "pars.a_st2", "pars.a_st3", "pars.a_y", "pars.beta1", "pars.curv_fact_colim", "pars.curv_fact_elec_trans", "pars.d", "pars.d_I", "pars.dmass_dN", "pars.eta", "pars.hmat", "pars.jmax_25", "pars.k_I", "pars.k_b", "pars.k_l", "pars.k_r", "pars.k_s", "pars.lma", "pars.nmass_b", "pars.nmass_l", "pars.nmass_r", "pars.nmass_s", "pars.omega", "pars.psi_crit", "pars.r_b", "pars.r_l", "pars.r_r", "pars.r_s", "pars.recruitment_decay", "pars.rho", "pars.root_b", "pars.root_c", "pars.root_depth_shape_eta", "pars.root_psi_crit", "pars.rooting_depth_max", "pars.stem_P50", "pars.stem_b", "pars.stem_c", "pars.theta", "pars.use_energy_balance", "pars.var_sapwood_volume_cost", "pars.vcmax_25"] + "value": ["birth_rate_x", "birth_rate_y", "control.GSS_tol_abs", "control.ci_abs_tol", "control.ci_niter", "control.fixed_time_step", "control.function_integration_rule", "control.node_density_in_birth_date", "control.node_gradient_direction", "control.node_gradient_eps", "control.node_gradient_richardson", "control.node_gradient_richardson_depth", "control.ode_a_dydt", "control.ode_a_y", "control.ode_step_size_initial", "control.ode_step_size_max", "control.ode_step_size_min", "control.ode_tol_abs", "control.ode_tol_rel", "control.offspring_production_iterations", "control.offspring_production_tol", "control.ppa_layer_optical_depth", "control.ppa_layer_smoothing", "control.save_RK45_cache", "control.schedule_eps", "control.schedule_nsteps", "control.schedule_verbose", "control.shading_model", "control.vulnerability_curve_ncontrol", "is_variable_birth_rate", "model_id", "pars.K_s", "pars.S_D", "pars.TF24_beta2", "pars.TF24_cost_scale", "pars.TF24_floor_lambda_o", "pars.a", "pars.a_b1", "pars.a_bio", "pars.a_cr1", "pars.a_d0", "pars.a_dG1", "pars.a_dG2", "pars.a_f1", "pars.a_f2", "pars.a_f3", "pars.a_l1", "pars.a_l2", "pars.a_p1", "pars.a_p2", "pars.a_r1", "pars.a_st1", "pars.a_st2", "pars.a_st3", "pars.a_y", "pars.beta1", "pars.curv_fact_colim", "pars.curv_fact_elec_trans", "pars.d", "pars.d_I", "pars.dmass_dN", "pars.eta", "pars.hmat", "pars.jmax_25", "pars.k_I", "pars.k_b", "pars.k_cr", "pars.k_l", "pars.k_r", "pars.k_s", "pars.lma", "pars.nmass_b", "pars.nmass_l", "pars.nmass_r", "pars.nmass_s", "pars.omega", "pars.psi_crit", "pars.r_b", "pars.r_cr", "pars.r_l", "pars.r_r", "pars.r_s", "pars.recruitment_decay", "pars.rho", "pars.root_b", "pars.root_c", "pars.root_depth_shape_eta", "pars.root_psi_crit", "pars.rooting_depth_max", "pars.stem_P50", "pars.stem_b", "pars.stem_c", "pars.theta", "pars.use_energy_balance", "pars.var_sapwood_volume_cost", "pars.vcmax_25"] } }, - "value": ["", "1", "0.001", "0.001", "1000", "0", "21", "FALSE", "-1", "9.9999999999999995e-07", "FALSE", "4", "0", "1", "9.9999999999999995e-07", "5", "9.9999999999999995e-07", "0.0001", "0.0001", "1000", "1e-08", "0.5", "0.29999999999999999", "FALSE", "0.02", "20", "FALSE", "", "100", "FALSE", "TF24f@v9.1", "1", "0.25", "1.5", "7.5", "0", "0.29999999999999999", "0.17000000000000001", "0.024500000000000001", "0.10000000000000001", "5.5", "20", "1", "50", "0.00011400000000000001", "5.4400000000000004", "0.30599999999999999", "151.177775377968", "0.204716166503633", "0.070000000000000007", "0.10000000000000001", "0.10000000000000001", "0.80000000000000004", "0.69999999999999996", "20000", "0.98999999999999999", "0.69999999999999996", "0.050000000000000003", "0.01", "0", "12", "16.595869100000002", "157.44", "0.5", "0.20000000000000001", "0.45658549999999998", "1", "0.20000000000000001", "0.1978791", "0.0033999999999999998", "0.012999999999999999", "0.0033500000000000001", "0.00198", "3.8000000000000002e-05", "7.0854928681315172", "13.197368421052632", "198.45451086041933", "217", "6.5986842105263159", "0", "608", "3.8982450000000002", "2.6801469999999998", "0.20000000000000001", "5.8702825428827037", "1.5", "1.8500000000000001", "2.5894373183081751", "1.0899851679673995", "0.00021417862497322766", "0", "1", "96"] + "value": ["", "1", "0.001", "0.001", "1000", "0", "21", "FALSE", "-1", "9.9999999999999995e-07", "FALSE", "4", "0", "1", "9.9999999999999995e-07", "5", "9.9999999999999995e-07", "0.0001", "0.0001", "1000", "1e-08", "0.5", "0.29999999999999999", "FALSE", "0.02", "20", "FALSE", "", "100", "FALSE", "TF24f@v9.1", "1", "0.25", "1.5", "7.5", "0", "0.29999999999999999", "0.17000000000000001", "0.024500000000000001", "0", "0.10000000000000001", "5.5", "20", "1", "50", "0.00011400000000000001", "5.4400000000000004", "0.30599999999999999", "151.177775377968", "0.204716166503633", "0.070000000000000007", "0.10000000000000001", "0.10000000000000001", "0.80000000000000004", "0.69999999999999996", "20000", "0.98999999999999999", "0.69999999999999996", "0.050000000000000003", "0.01", "0", "12", "16.595869100000002", "157.44", "0.5", "0.20000000000000001", "0.20000000000000001", "0.45658549999999998", "1", "0.20000000000000001", "0.1978791", "0.0033999999999999998", "0.012999999999999999", "0.0033500000000000001", "0.00198", "3.8000000000000002e-05", "7.0854928681315172", "13.197368421052632", "6.5986842105263159", "198.45451086041933", "217", "6.5986842105263159", "0", "608", "3.8982450000000002", "2.6801469999999998", "0.20000000000000001", "5.8702825428827037", "1.5", "1.8500000000000001", "2.5894373183081751", "1.0899851679673995", "0.00021417862497322766", "0", "1", "96"] } diff --git a/tests/testthat/test-coarse-roots.R b/tests/testthat/test-coarse-roots.R new file mode 100644 index 000000000..d9cf41108 --- /dev/null +++ b/tests/testthat/test-coarse-roots.R @@ -0,0 +1,181 @@ +## Coarse (structural) root pool -- issue #349. +## +## FF16 and TF24 gained a second root pool, m_cr = a_cr1 * m_s: coarse roots as +## the below-ground continuation of the sapwood cylinder, so they inherit its +## height scaling and the structural root:shoot ratio is size-invariant. Fine +## roots (a_r1 * A_l) are unchanged. +## +## The pool ships OFF (a_cr1 = 0), which splits these tests in two. Half of them +## check it really is inert at the defaults -- the reason FF16 keeps agreeing +## with the published 2012 model in test-strategy-ff16-reference-comparison.R +## and no scientific version moved. The other half switch it on, because a +## default-off feature only ever run at its default is not tested at all. + +eta_c_of <- function(p) 1 - 2 / (1 + p$eta) + 1 / (1 + 2 * p$eta) + +## Published FF16 allometric derivatives, written out here rather than read back +## from the strategy, so the C++ is compared against the documented formulas +## (inst/docs/FF16/FF16-eqns.csv) and not against itself. +dmass_sapwood_darea_leaf_r <- function(p, area_leaf) { + p$rho * eta_c_of(p) * p$a_l1 * p$theta * (p$a_l2 + 1) * area_leaf^p$a_l2 +} + +area_leaf_r <- function(p, height) (height / p$a_l1)^(1 / p$a_l2) + +## Sum of d(component mass)/d(area_leaf) over the live pools -- the denominator +## of darea_leaf_dmass_live, i.e. the mass that must be built alongside each new +## unit of leaf area. +dmass_live_darea_leaf_r <- function(p, height) { + dms <- dmass_sapwood_darea_leaf_r(p, area_leaf_r(p, height)) + p$lma + dms + p$a_b1 * dms + p$a_r1 + p$a_cr1 * dms +} + +allometry_of <- function(type, s, height) { + z <- rep(0, length(height)) + switch(type, + FF16 = FF16_strategy_expand_allometry(s, height, z, z), + TF24 = TF24_strategy_expand_allometry(s, height, z, z)) +} + +heights <- seq(0.5, 20, length.out = 40) + +test_that("the coarse-root pool is absent at the shipped defaults", { + for (type in c("FF16", "TF24")) { + s <- get(paste0(type, "_Strategy"))() + expect_identical(s$pars$a_cr1, 0, info = type) + + al <- allometry_of(type, s, heights) + expect_true(all(al$mass_coarse_root == 0), info = type) + ## Live and total mass therefore still read exactly as they did before the + ## pool existed -- appending a term worth 0 changes no bits. + expect_identical(al$mass_live, + al$mass_leaf + al$mass_sapwood + al$mass_bark + al$mass_root, + info = type) + } +}) + +test_that("coarse-root mass is the sapwood fraction, and the mass budget closes", { + a_cr1 <- 0.25 + for (type in c("FF16", "TF24")) { + s <- get(paste0(type, "_Strategy"))() + s$pars$a_cr1 <- a_cr1 + ## Non-zero heartwood, so mass_total and mass_above_ground are not just + ## restatements of mass_live. + mh <- seq(0, 5, length.out = length(heights)) + ah <- seq(0, 0.01, length.out = length(heights)) + al <- switch(type, + FF16 = FF16_strategy_expand_allometry(s, heights, ah, mh), + TF24 = TF24_strategy_expand_allometry(s, heights, ah, mh)) + + expect_equal(al$mass_coarse_root, a_cr1 * al$mass_sapwood, + tolerance = 1e-14, info = type) + ## Coarse roots are below ground: they belong to the total but not to the + ## above-ground mass, and the two root pools are the whole difference. + expect_equal(al$mass_total, + al$mass_above_ground + al$mass_root + al$mass_coarse_root, + tolerance = 1e-12, info = type) + ## Size-invariant structural root:shoot -- the property that motivated + ## scaling on sapwood rather than on leaf area. + expect_equal(al$mass_coarse_root / al$mass_sapwood, + rep(a_cr1, length(heights)), tolerance = 1e-14, info = type) + } +}) + +test_that("coarse roots slow height growth by exactly the allocation ratio", { + ## Setting r_cr = k_cr = 0 makes the pool free to maintain, so net production + ## is untouched and the ONLY route left is allocation: each new unit of leaf + ## area now has to be built with a_cr1 * dm_s/dA_l of extra coarse root. The + ## predicted slowdown is therefore the ratio of the two live-mass denominators, + ## which is a closed form the strategy never computes. + a_cr1 <- 0.3 + h0 <- 8 + + dh_dt <- function(a_cr1, r_cr, k_cr) { + s <- FF16_Strategy() + s$pars$a_cr1 <- a_cr1 + s$pars$r_cr <- r_cr + s$pars$k_cr <- k_cr + env <- Environment("FF16") + env$set_fixed_environment(1.0, height_max = 150) + pl <- FF16_Individual(s) + pl$set_state("height", h0) + pl$compute_rates(env) + list(dh = pl$rate("height"), P = pl$aux("net_mass_production_dt"), pars = s$pars) + } + + off <- dh_dt(0, 0, 0) + on <- dh_dt(a_cr1, 0, 0) + + ## Free maintenance => identical carbon income. + expect_identical(on$P, off$P) + + predicted <- dmass_live_darea_leaf_r(off$pars, h0) / + dmass_live_darea_leaf_r(on$pars, h0) + expect_equal(on$dh / off$dh, predicted, tolerance = 1e-12) + expect_lt(on$dh, off$dh) + + ## The ratio alone would survive an error common to both denominators, so + ## reconstruct dh/dt outright from the published chain + ## dH/dt = dH/dA_l * P * (1 - r(H)) / (dM_live/dA_l) + ## and check the absolute value with the pool switched on. + dheight_darea_leaf_r <- function(p, height) { + p$a_l1 * p$a_l2 * area_leaf_r(p, height)^(p$a_l2 - 1) + } + frac_growth_r <- function(p, height) { + 1 - p$a_f1 / (1 + exp(p$a_f2 * (1 - height / p$hmat))) + } + expect_equal(on$dh, + dheight_darea_leaf_r(on$pars, h0) * on$P * + frac_growth_r(on$pars, h0) / + dmass_live_darea_leaf_r(on$pars, h0), + tolerance = 1e-12) +}) + +test_that("coarse roots also cost carbon to maintain", { + ## At the shipped woody rates (r_cr = r_s, k_cr = k_s) the pool costs on both + ## sides: respiration and turnover cut net production as well. So switching it + ## on must slow growth strictly more than the allocation effect alone. + h0 <- 8 + rates <- function(a_cr1, free) { + s <- FF16_Strategy() + s$pars$a_cr1 <- a_cr1 + if (free) { s$pars$r_cr <- 0; s$pars$k_cr <- 0 } + env <- Environment("FF16") + env$set_fixed_environment(1.0, height_max = 150) + pl <- FF16_Individual(s) + pl$set_state("height", h0) + pl$compute_rates(env) + c(dh = pl$rate("height"), P = pl$aux("net_mass_production_dt")) + } + + off <- rates(0.0, FALSE) + alloc <- rates(0.3, TRUE) + alloc_mnt <- rates(0.3, FALSE) + + expect_lt(alloc_mnt[["P"]], off[["P"]]) # maintenance eats production + expect_lt(alloc_mnt[["dh"]], alloc[["dh"]]) # ... on top of the allocation cost + expect_lt(alloc[["dh"]], off[["dh"]]) +}) + +test_that("coarse roots stay out of TF24's water uptake", { + ## TF24 distributes FINE-root mass over soil depth to drive uptake, and reads + ## a_r1 directly to do it. Coarse roots are structure, not absorbing surface, + ## so at a fixed size every water-side quantity must be bit-identical whether + ## the pool is on or off. + water_aux <- function(a_cr1) { + s <- TF24_Strategy() + s$pars$a_cr1 <- a_cr1 + env <- Environment("TF24") + env$set_fixed_environment(1.0, height_max = 150) + env$set_soil_water_state(rep(0.4, env$get_soil_number_of_depths())) + env$time <- 5 + ind <- Individual("TF24", "TF24_Env")(s) + ind$set_state("height", 8) + ind$compute_rates(env) + vapply(c("root_mass", "transpiration", "E_up_", "opt_psi_stem", + "opt_root_psi", "assimilation"), + function(nm) ind$aux(nm), numeric(1)) + } + + expect_identical(water_aux(0.3), water_aux(0.0)) +}) diff --git a/tests/testthat/test-expand-state.R b/tests/testthat/test-expand-state.R index 2fb293208..0f67273a7 100644 --- a/tests/testthat/test-expand-state.R +++ b/tests/testthat/test-expand-state.R @@ -15,18 +15,21 @@ allometry_oracle <- function(s, height, area_heartwood, mass_heartwood) { area_stem <- area_bark + area_sapwood + area_heartwood diameter_stem <- sqrt(4 * area_stem / pi) mass_root <- p$a_r1 * area_leaf - mass_live <- mass_leaf + mass_sapwood + mass_bark + mass_root - mass_total <- mass_leaf + mass_bark + mass_sapwood + mass_heartwood + mass_root + mass_coarse_root <- p$a_cr1 * mass_sapwood + mass_live <- mass_leaf + mass_sapwood + mass_bark + mass_root + mass_coarse_root + mass_total <- mass_leaf + mass_bark + mass_sapwood + mass_heartwood + mass_root + mass_coarse_root mass_above_ground <- mass_leaf + mass_bark + mass_sapwood + mass_heartwood list(area_leaf = area_leaf, mass_leaf = mass_leaf, area_sapwood = area_sapwood, mass_sapwood = mass_sapwood, area_bark = area_bark, mass_bark = mass_bark, area_stem = area_stem, diameter_stem = diameter_stem, mass_root = mass_root, + mass_coarse_root = mass_coarse_root, mass_live = mass_live, mass_total = mass_total, mass_above_ground = mass_above_ground) } -check_allometry <- function(strategy_fn, cpp_fn) { +check_allometry <- function(strategy_fn, cpp_fn, a_cr1 = 0) { s <- strategy_fn() + s$pars$a_cr1 <- a_cr1 height <- seq(0.5, 20, length.out = 50) # vary heartwood so area_stem / mass_total / mass_above_ground are exercised area_heartwood <- seq(0, 0.01, length.out = 50) @@ -49,6 +52,18 @@ test_that("TF24 C++ allometry matches the historical R formulas", { check_allometry(TF24_Strategy, TF24_strategy_expand_allometry) }) +## The oracle above carries the coarse-root pool (#349), but at the shipped +## default a_cr1 = 0 every coarse-root term is zero, so the checks above would +## pass even if C++ ignored a_cr1 entirely. Re-run them with the pool switched +## on so the new terms are actually compared. +test_that("FF16 C++ allometry matches the R formulas with coarse roots on", { + check_allometry(FF16_Strategy, FF16_strategy_expand_allometry, a_cr1 = 0.25) +}) + +test_that("TF24 C++ allometry matches the R formulas with coarse roots on", { + check_allometry(TF24_Strategy, TF24_strategy_expand_allometry, a_cr1 = 0.25) +}) + test_that("FF16_expand_state adds the expected derived columns", { p <- scm_base_parameters("FF16") p <- add_strategies(p, trait_matrix(0.0825, "lma"), hyperpar = FF16_hyperpar, birth_rate = list(20)) @@ -57,7 +72,8 @@ test_that("FF16_expand_state adds the expected derived columns", { derived <- c("area_leaf", "mass_leaf", "area_sapwood", "mass_sapwood", "area_bark", "mass_bark", "area_stem", "diameter_stem", - "mass_root", "mass_live", "mass_total", "mass_above_ground") + "mass_root", "mass_coarse_root", "mass_live", "mass_total", + "mass_above_ground") expect_true(all(derived %in% names(es$species))) # Re-derive with the oracle from the (already present) state columns and the diff --git a/tests/testthat/test-ff16-ad-kernel.R b/tests/testthat/test-ff16-ad-kernel.R index 55075cab3..e72d65421 100644 --- a/tests/testthat/test-ff16-ad-kernel.R +++ b/tests/testthat/test-ff16-ad-kernel.R @@ -47,6 +47,7 @@ compile_ff16_ad_kernel <- function() { p.lma=v[0];p.rho=v[1];p.theta=v[2];p.a_b1=v[3];p.a_r1=v[4];p.eta_c=v[5]; p.a_p1=v[6];p.a_p2=v[7];p.r_l=v[8];p.r_s=v[9];p.r_b=v[10];p.r_r=v[11]; p.k_l=v[12];p.k_b=v[13];p.k_s=v[14];p.k_r=v[15];p.a_bio=v[16];p.a_y=v[17]; + p.a_cr1=v[18];p.r_cr=v[19];p.k_cr=v[20]; return p; } static plant::FF16ProdPars pod_d(const std::vector& v) { @@ -54,6 +55,7 @@ compile_ff16_ad_kernel <- function() { p.lma=v[0];p.rho=v[1];p.theta=v[2];p.a_b1=v[3];p.a_r1=v[4];p.eta_c=v[5]; p.a_p1=v[6];p.a_p2=v[7];p.r_l=v[8];p.r_s=v[9];p.r_b=v[10];p.r_r=v[11]; p.k_l=v[12];p.k_b=v[13];p.k_s=v[14];p.k_r=v[15];p.a_bio=v[16];p.a_y=v[17]; + p.a_cr1=v[18];p.r_cr=v[19];p.k_cr=v[20]; return p; } @@ -94,16 +96,22 @@ testthat::test_that("FF16 net-production AD gradient matches finite differences" "Skipping FF16 AD kernel in pkgload load_all sessions.") compile_ff16_ad_kernel() + ## Trailing three are the coarse-root pool (#349): a_cr1, r_cr, k_cr. a_cr1 is + ## deliberately non-zero here so the new terms are actually differentiated -- + ## at the shipped default of 0 the whole pool drops out of the gradient. v <- c(0.1978791, 608, 0.0002141786, 0.17, 0.07, 0.5805, 151.177, 0.204, - 0.01979, 0.0859, 0.04, 0.2086, 0.4565, 0.2, 0.0, 1.0, 0.0245, 0.7) + 0.01979, 0.0859, 0.04, 0.2086, 0.4565, 0.2, 0.0, 1.0, 0.0245, 0.7, + 0.2, 0.0859, 0.2) height <- 5; area_leaf <- 0.3; light_E <- 0.78 g <- ff16_netprod_grad(v, height, area_leaf, light_E) inp <- c(v, height, area_leaf, light_E) + np <- length(v) fd <- vapply(seq_along(inp), function(i) { h <- 1e-6 * max(1, abs(inp[i])) up <- inp; dn <- inp; up[i] <- up[i] + h; dn[i] <- dn[i] - h - f <- function(z) ff16_netprod_value(z[1:18], z[19], z[20], z[21]) + f <- function(z) ff16_netprod_value(z[seq_len(np)], z[np + 1], z[np + 2], + z[np + 3]) (f(up) - f(dn)) / (2 * h) }, numeric(1)) diff --git a/tests/testthat/test-ff16-resident-coupling-ad.R b/tests/testthat/test-ff16-resident-coupling-ad.R index 32d108533..c30b22906 100644 --- a/tests/testthat/test-ff16-resident-coupling-ad.R +++ b/tests/testthat/test-ff16-resident-coupling-ad.R @@ -65,6 +65,7 @@ compile_ff16_resident_ad <- function() { p.lma=v[0];p.rho=v[1];p.theta=v[2];p.a_b1=v[3];p.a_r1=v[4];p.eta_c=v[5]; p.a_p1=v[6];p.a_p2=v[7];p.r_l=v[8];p.r_s=v[9];p.r_b=v[10];p.r_r=v[11]; p.k_l=v[12];p.k_b=v[13];p.k_s=v[14];p.k_r=v[15];p.a_bio=v[16];p.a_y=v[17]; + p.a_cr1=v[18];p.r_cr=v[19];p.k_cr=v[20]; return p; } @@ -111,8 +112,11 @@ testthat::test_that("FF16 resident self-shading coupling differentiates end-to-e wq <- simp * dz / 3 * q a_l1 <- 0.306; a_l2 <- 0.75; k_I <- 0.5 + ## Trailing three are the coarse-root pool (#349): a_cr1, r_cr, k_cr, with a + ## non-zero a_cr1 so the a_l1 gradient runs through the new terms too. v <- c(0.1978791, 608, 0.0002141786, 0.17, 0.07, 0.5805, 151.177, 0.204, - 0.01979, 0.0859, 0.04, 0.2086, 0.4565, 0.2, 0.0, 1.0, 0.0245, 0.7) + 0.01979, 0.0859, 0.04, 0.2086, 0.4565, 0.2, 0.0, 1.0, 0.0245, 0.7, + 0.2, 0.0859, 0.2) g <- resident_grad(a_l1, k_I, v, a_l2, H, density, xk, qfrac_knot, z, wq) f <- function(al1, ki) resident_value(al1, ki, v, a_l2, H, density, xk, diff --git a/tests/testthat/test-strategy-ff16-reference-comparison.R b/tests/testthat/test-strategy-ff16-reference-comparison.R index bffb7f69f..9761ca548 100644 --- a/tests/testthat/test-strategy-ff16-reference-comparison.R +++ b/tests/testthat/test-strategy-ff16-reference-comparison.R @@ -14,8 +14,11 @@ test_that("FF16_Strategy parameters agree with reference model", { v <- setdiff(names(cmp_pars), c("n_area")) expect_true(all(v %in% names(sp))) - ## And v.v., except for a few additions present in pars but not the reference: - extra <- c("S_D", "recruitment_decay") + ## And v.v., except for a few additions present in pars but not the reference. + ## The coarse-root trio (#349) is absent from the 2012 reference model by + ## construction: at the shipped default a_cr1 = 0 the pool contributes nothing, + ## so the two implementations still agree on every growth quantity below. + extra <- c("S_D", "recruitment_decay", "a_cr1", "r_cr", "k_cr") common <- setdiff(names(sp), extra) expect_true(all(extra %in% names(sp))) expect_true(all(common %in% names(cmp_pars))) diff --git a/tests/testthat/test-strategy-ff16.R b/tests/testthat/test-strategy-ff16.R index d005295f2..55747fd78 100644 --- a/tests/testthat/test-strategy-ff16.R +++ b/tests/testthat/test-strategy-ff16.R @@ -7,8 +7,10 @@ test_that("Defaults", { a_y = 0.7, a_l1 = 5.44, a_r1 = 0.07, + a_cr1 = 0, a_b1 = 0.17, r_b = 8024 / 608, + r_cr = 4012/608, r_l = 39.27 / 0.1978791, r_r = 217, r_s = 4012/608, @@ -26,6 +28,7 @@ test_that("Defaults", { hmat = 16.5958691, k_b = 0.2, k_l = 0.4565855, + k_cr = 0.2, k_r = 1, k_s = 0.2, lma = 0.1978791, diff --git a/tests/testthat/test-strategy-tf24.R b/tests/testthat/test-strategy-tf24.R index f321f1468..5b4428602 100644 --- a/tests/testthat/test-strategy-tf24.R +++ b/tests/testthat/test-strategy-tf24.R @@ -8,8 +8,10 @@ test_that("Defaults", { a_y = 0.7, a_l1 = 5.44, a_r1 = 0.07, + a_cr1 = 0, a_b1 = 0.17, r_b = 8024 / 608, + r_cr = 4012/608, r_l = 39.27 / 0.1978791, r_r = 217, r_s = 4012/608, @@ -30,6 +32,7 @@ test_that("Defaults", { hmat = 16.5958691, k_b = 0.2, k_l = 0.4565855, + k_cr = 0.2, k_r = 1, k_s = 0.2, lma = 0.1978791, diff --git a/tests/testthat/test-strategy-tf24f.R b/tests/testthat/test-strategy-tf24f.R index 2c22c5cb9..c81cdc56c 100644 --- a/tests/testthat/test-strategy-tf24f.R +++ b/tests/testthat/test-strategy-tf24f.R @@ -9,8 +9,10 @@ test_that("Defaults", { a_y = 0.7, a_l1 = 5.44, a_r1 = 0.07, + a_cr1 = 0, a_b1 = 0.17, r_b = 8024 / 608, + r_cr = 4012/608, r_l = 39.27 / 0.1978791, r_r = 217, r_s = 4012/608, @@ -31,6 +33,7 @@ test_that("Defaults", { hmat = 16.5958691, k_b = 0.2, k_l = 0.4565855, + k_cr = 0.2, k_r = 1, k_s = 0.2, lma = 0.1978791,