Calculation of annual fAPAR and LAI values from a fitted PModel and precipitation data

[davidorme]

This is the first step in implementing #347 and the description below is taken from the draft implementation document in the meta issue.

Theory

The predicted fAPARmax / LAImax mainly has two roles:

1. is to calculate a parameter m which is the ratio of steady-state LAI to steady-state GPP;
2. is to limit the prediction of LAI timeseries, predicted LAI ≤ LAImax.

The equations to calculate fAPARmax (Eqn 1) and LAImax (Eqn 2) are:

(Eqn. 1) $f{APAR{max}}= min\left[ (1 – z/(kA_0)), [c_a (1–\chi)/1.6 D] [f_0 P/A0]\right]$ (Eqn. 2) $LAI{max}= –(1/k) ln (1 –f{APAR{opt}})$

Of those terms, the following are new constants across the calculations. There aren't too many, so I think these are parameters to the calculation but we might want to bundle them into a new constants class (MaxFAPARConst?)

• $z$ accounts for the costs of building and maintaining leaves and the total below-ground allocation required to support the nutrient demand of those leaves, and the globally fitted values for z was 12.227 mol m–2 year–1.
• $k$ is the light extinction coefficient (set at 0.5).
• $f_0$ is the ratio of annual total transpiration of annual total precipitation, which is an empirical function (Eqn 3) of the climatic Aridity Index (AI).

(Eqn 3) $f_0=0.65 e^{-b ln^2 \left(\frac{AI}{1.9}\right) }$

where 0.65 is the maximum value, 1.9 is the AI at which this maximum occurs, and b = 0.604169. For AI , see below.

The following need to be be extracted from a P Model

• $A_0$ is the annual sum of potential GPP. Potential GPP would be achieved if fAPAR = 1. (mol m-2 year-1).
• $c_a$ is the ambient CO2 partial pressure (Pa).
• $\chi$ is the annual mean ratio of leaf-internal CO2 partial pressure (Pa) to ca during the growing season (>0℃).
• $D$ is the annual mean vapour pressure deficit (VPD, Pa) during the growing season (>0℃).

And lastly, precipitation data and aridity data is needed:

• P is the annual total precipitation (mol m–2 year–1), *AI is a climatological estimate of local aridity index, typically calculated as total PET over total precipitation for an appropriate period, typically at least 20 years.

An Example to predict fAPARmax:

z is 12.227 mol m–2 year–1; 
k is 0.5; 
A0 is 200 mol m–2 year–1; 
ca is 40 Pa;
χ is 0.8; 
D is 3500 Pa; 
P is 51888 mol m–2 year–1; 
f0 is 0.62.
fAPAR_max=min⁡{1 –12.227/(0.5×200),[37×(1–0.8)/(1.6×3500)][0.62×51888/200]}  
                  = min {0.87, 0.21}
                  = 0.21

Implementation

This is a draft workflow for implementing the calculation

• [ ] A user fits a P Model - this could be the standard or subdaily model. The model must cover at least one year and should be fitted with a $fAPAR$ fixed at 1 to give predictions of potential GPP ($A_0$ above). This could be fitted to a single site or to arrays for multiple sites.
• [ ] The user needs to provide a sequence of datetimes for each observation, which have to map onto one of the axes of the input data. For now, we follow the model used in the SubdailyPModel and SPLASH of insisting this is the first axis of the arrays.
• [ ] Get the growing seasons from the datetimes and temperatures (see #352).
• [ ] Using the PModel, the datetimes and the growing season, we can then extract annual values for annual mean $c_a, \chi, D$.
• [ ] We then need to take user provided precipitation to get annual precipitation.

So, I think the signature will look like this:

class MaximumFAPAR:

    def __init__(
        self, 
        pmodel: PModel, 
        datetimes: NPArray[np.datetime64],
        growing_season: NPArray[np.bool_],
        precipitation: NPArray[np.datetime64],
        aridity_index: float,
        additional_constants: ...
    ):
        ...

The requirement for precipitation and aridity index ties this quite closely to the SPLASH module. The demonstration implementation should probably use SPLASH to set up the scenario.

[zhouboya]

For the parameter z, insteady of directly giving a global fitted value (12.227 mol m-2 year-1), Shirely suggested it is better to allow the users to fit this parameter values based on their own satellite/climate data.

z accounts for the costs of building and maintaining leaves and the total below-ground allocation required to support the nutrient demand of those leaves and should have a globally fitted value for simplicity. Considering the fitted z value could vary with inputted satellite data, thus, it is better to fit z value based on the satellite data users have.

Methods to fit z values are:

In equation 1, fAPARmax is calculated as the minimum values of energy-limited fAPAR (fAPARc, Eqn5) and water-limited fAPAR (fAPARw, Eqn 4).

(Eqn. 4) fAPARw= [𝑐𝑎 (1–𝜒)/1.6 𝐷] [𝑓0 𝑃/𝐴0] (Eqn. 5) fAPARc= 1 – 𝑧/(𝑘𝐴0)

The minimum of energy-limited fAPAR (fAPARc, Eqn 5) and water-limited fAPAR (fAPARw, Eqn4) was smoothly approximated using the log-sum-exp function:

(Eqn.6) 𝑚𝑖𝑛{𝑓𝐴𝑃𝐴𝑅𝑐,𝑓𝐴𝑃𝐴𝑅𝑤}≈−1𝐾𝑙𝑜𝑔[exp(−𝐾.𝑓𝐴𝑃𝐴𝑅𝑐)+exp(−𝐾.𝑓𝐴𝑃𝐴𝑅𝑤)]

where K (≥ 1) is a constant. The larger the value for K, the closer the approximation is to minimum function (here we adopted K=10). With fAPARw and fAPARc defined by equations (4) and (5) respectively, parameter z can be fitted using non-linear regression (nls function in R), with user's own satellite data as inputs.

[zhouboya]

For the inputted data (P) of Eqn 1: P, which is the annual total precipitation (mol m–2 year–1). One thing needs to be noticed is that, here the units of P is mol m-2 year-1. But for most precipitation data, their given units are mm. For users, they may often overlook this point or make mistakes easily. So it is better to include a unit converter in the package to convert mm to mol, for the convenience of user.