Add an optional coarse-root pool to FF16 and TF24

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

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

inst/RcppR6_classes.yml

@@ -1150,17 +1150,20 @@ 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
     - k_l: double
     - k_b: double
     - k_s: double
     - k_r: double
+    - k_cr: double
     - a_p1: double
     - a_p2: double
     - a_f3: double
@@ -1319,17 +1322,20 @@ 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
     - k_l: double
     - k_b: double
     - k_s: double
     - k_r: double
+    - k_cr: double
     - a_p1: double
     - a_p2: double
     - a_f3: double

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}$"
 

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

inst/include/plant/RcppR6_post.hpp

@@ -1105,17 +1105,20 @@ 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);
   ret["k_l"] = Rcpp::wrap(x.k_l);
   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<double >(xl["a_l2"]);
   // ret.a_r1 = Rcpp::as<decltype(reta_r1) >(xl["a_r1"]);
   ret.a_r1 = Rcpp::as<double >(xl["a_r1"]);
+  // ret.a_cr1 = Rcpp::as<decltype(reta_cr1) >(xl["a_cr1"]);
+  ret.a_cr1 = Rcpp::as<double >(xl["a_cr1"]);
   // ret.a_b1 = Rcpp::as<decltype(reta_b1) >(xl["a_b1"]);
   ret.a_b1 = Rcpp::as<double >(xl["a_b1"]);
   // ret.r_s = Rcpp::as<decltype(retr_s) >(xl["r_s"]);
@@ -1166,6 +1171,8 @@ template <> inline plant::FF16_Pars as(SEXP x) {
   ret.r_b = Rcpp::as<double >(xl["r_b"]);
   // ret.r_r = Rcpp::as<decltype(retr_r) >(xl["r_r"]);
   ret.r_r = Rcpp::as<double >(xl["r_r"]);
+  // ret.r_cr = Rcpp::as<decltype(retr_cr) >(xl["r_cr"]);
+  ret.r_cr = Rcpp::as<double >(xl["r_cr"]);
   // ret.r_l = Rcpp::as<decltype(retr_l) >(xl["r_l"]);
   ret.r_l = Rcpp::as<double >(xl["r_l"]);
   // ret.a_y = Rcpp::as<decltype(reta_y) >(xl["a_y"]);
@@ -1180,6 +1187,8 @@ template <> inline plant::FF16_Pars as(SEXP x) {
   ret.k_s = Rcpp::as<double >(xl["k_s"]);
   // ret.k_r = Rcpp::as<decltype(retk_r) >(xl["k_r"]);
   ret.k_r = Rcpp::as<double >(xl["k_r"]);
+  // ret.k_cr = Rcpp::as<decltype(retk_cr) >(xl["k_cr"]);
+  ret.k_cr = Rcpp::as<double >(xl["k_cr"]);
   // ret.a_p1 = Rcpp::as<decltype(reta_p1) >(xl["a_p1"]);
   ret.a_p1 = Rcpp::as<double >(xl["a_p1"]);
   // ret.a_p2 = Rcpp::as<decltype(reta_p2) >(xl["a_p2"]);
@@ -1349,17 +1358,20 @@ 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);
   ret["k_l"] = Rcpp::wrap(x.k_l);
   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<double >(xl["a_l2"]);
   // ret.a_r1 = Rcpp::as<decltype(reta_r1) >(xl["a_r1"]);
   ret.a_r1 = Rcpp::as<double >(xl["a_r1"]);
+  // ret.a_cr1 = Rcpp::as<decltype(reta_cr1) >(xl["a_cr1"]);
+  ret.a_cr1 = Rcpp::as<double >(xl["a_cr1"]);
   // ret.a_b1 = Rcpp::as<decltype(reta_b1) >(xl["a_b1"]);
   ret.a_b1 = Rcpp::as<double >(xl["a_b1"]);
   // ret.r_s = Rcpp::as<decltype(retr_s) >(xl["r_s"]);
@@ -1440,6 +1454,8 @@ template <> inline plant::TF24_Pars as(SEXP x) {
   ret.r_b = Rcpp::as<double >(xl["r_b"]);
   // ret.r_r = Rcpp::as<decltype(retr_r) >(xl["r_r"]);
   ret.r_r = Rcpp::as<double >(xl["r_r"]);
+  // ret.r_cr = Rcpp::as<decltype(retr_cr) >(xl["r_cr"]);
+  ret.r_cr = Rcpp::as<double >(xl["r_cr"]);
   // ret.r_l = Rcpp::as<decltype(retr_l) >(xl["r_l"]);
   ret.r_l = Rcpp::as<double >(xl["r_l"]);
   // ret.a_y = Rcpp::as<decltype(reta_y) >(xl["a_y"]);
@@ -1454,6 +1470,8 @@ template <> inline plant::TF24_Pars as(SEXP x) {
   ret.k_s = Rcpp::as<double >(xl["k_s"]);
   // ret.k_r = Rcpp::as<decltype(retk_r) >(xl["k_r"]);
   ret.k_r = Rcpp::as<double >(xl["k_r"]);
+  // ret.k_cr = Rcpp::as<decltype(retk_cr) >(xl["k_cr"]);
+  ret.k_cr = Rcpp::as<double >(xl["k_cr"]);
   // ret.a_p1 = Rcpp::as<decltype(reta_p1) >(xl["a_p1"]);
   ret.a_p1 = Rcpp::as<double >(xl["a_p1"]);
   // ret.a_p2 = Rcpp::as<decltype(reta_p2) >(xl["a_p2"]);

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 <typename S>
 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 <typename S>
 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 <typename S>
 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<S>& 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);
 }
 

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<FF16_Environment> {
   // [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<FF16_Environment> {
 
   // [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<FF16_Environment> {
   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<FF16_Environment> {
   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,