Documentation fixes from Giulia Mengoli and updated refs

docs/source/refs.bib

@@ -732,3 +732,109 @@ @article{mengoli:2023a
 	pages = {1--19},
 	file = {Full Text PDF:/Users/dorme/Zotero/storage/JVIB9NS3/Mengoli et al. - 2023 - A global function of climatic aridity accounts for.pdf:application/pdf},
 }
+
+@techreport{allen:1998a,
+	address = {Rome},
+	title = {Crop evapotranspiration - {Guidelines} for computing crop water requirements},
+	number = {56},
+	institution = {FAO},
+	author = {Allen, Richard G. and Pereira, Luis S and Raes, Dirk and Smith, Martin},
+	year = {1998},
+	file = {Allen et al. - 1998 - Crop evapotranspiration - Guidelines for computing.pdf:/Users/dorme/Zotero/storage/68AVA64R/Allen et al. - 1998 - Crop evapotranspiration - Guidelines for computing.pdf:application/pdf},
+}
+
+@article{tsilingiris:2008a,
+	title = {Thermophysical and transport properties of humid air at temperature range between 0 and 100°{C}},
+	volume = {49},
+	issn = {0196-8904},
+	url = {https://www.sciencedirect.com/science/article/pii/S0196890407003329},
+	doi = {10.1016/j.enconman.2007.09.015},
+	abstract = {The aim of the present investigation is evaluation of the thermophysical and transport properties of moist air as a function of mixture temperature with relative humidity as a parameter, ranging between dry air and saturation conditions. Based on a literature review of the most widely available analytical procedures and methods, a number of developed correlations are presented, which are employed with recent gas mixture component properties as input parameters, to derive the temperature and humidity dependence of mixture density, viscosity, specific heat capacity, thermal conductivity, thermal diffusivity and Prandtl number under conditions corresponding to the total barometric pressure of 101.3kPa. The derived results at an accuracy level suitable for engineering calculations were plotted and compared with adequate accuracy with existing results from previous analytical calculations and measured data from earlier experimental investigations. The saturated mixture properties were also appropriately fitted, and the fitting expressions suitable for computer calculations are also presented.},
+	language = {en},
+	number = {5},
+	urldate = {2023-07-04},
+	journal = {Energy Conversion and Management},
+	author = {Tsilingiris, P. T.},
+	month = may,
+	year = {2008},
+	keywords = {Density, Prandtl number, Specific heat capacity, Thermal conductivity, Thermal diffusivity, Thermophysical properties, Viscosity},
+	pages = {1098--1110},
+	file = {ScienceDirect Full Text PDF:/Users/dorme/Zotero/storage/EZSDLB2V/Tsilingiris - 2008 - Thermophysical and transport properties of humid a.pdf:application/pdf;ScienceDirect Snapshot:/Users/dorme/Zotero/storage/IZWWYPC8/S0196890407003329.html:text/html},
+}
+
+@article{henderson-sellers:1984a,
+	title = {A new formula for latent heat of vaporization of water as a function of temperature},
+	volume = {110},
+	copyright = {Copyright © 1984 Royal Meteorological Society},
+	issn = {1477-870X},
+	url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/qj.49711046626},
+	doi = {10.1002/qj.49711046626},
+	abstract = {Existing formulae and approximations for the latent heat of vaporization of water, Lv, are reviewed. Using an analytical approximation to the saturated vapour pressure as a function of temperature, a new, temperature-dependent function for Lv is derived.},
+	language = {en},
+	number = {466},
+	urldate = {2023-07-04},
+	journal = {Quarterly Journal of the Royal Meteorological Society},
+	author = {Henderson-Sellers, B.},
+	year = {1984},
+	note = {\_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/qj.49711046626},
+	pages = {1186--1190},
+	file = {Snapshot:/Users/dorme/Zotero/storage/WJAU4R6F/qj.html:text/html},
+}
+
+@article{chen:2008a,
+	title = {The equation of state of pure water determined from sound speeds},
+	volume = {66},
+	issn = {0021-9606},
+	url = {https://doi.org/10.1063/1.434179},
+	doi = {10.1063/1.434179},
+	abstract = {The equation of state of water valid over the range 0–100 °C and 0–1000 bar has been determined from the high pressure sound velocities of Wilson, which were reanalyzed by Chen and Millero. The equation of state has a maximum error of ±0.01 bar−1 in isothermal compressibility and is in the form of a secant bulk modulus: K=V0P/(V0−V) =K0+AP+BP2, where K, K0, and V, V0 are the secant bulk moduli and specific volumes at applied pressures P and 0 (1 atm), respectively; A and B are temperature dependent parameters. The good agreement (to within 20×10−6 cm3 g−1) of specific volumes calculated using the above equation with those obtained from other modifications of the Wilson sound velocity data demonstrates the reliability of the sound velocity method for determining equations of state.},
+	number = {5},
+	urldate = {2023-07-04},
+	journal = {The Journal of Chemical Physics},
+	author = {Chen, Chen‐Tung and Fine, Rana A. and Millero, Frank J.},
+	month = aug,
+	year = {2008},
+	pages = {2142--2144},
+	file = {Full Text PDF:/Users/dorme/Zotero/storage/2WBI492T/Chen et al. - 2008 - The equation of state of pure water determined fro.pdf:application/pdf;Snapshot:/Users/dorme/Zotero/storage/IRHA5IDN/The-equation-of-state-of-pure-water-determined.html:text/html},
+}
+
+@article{hengl:2017a,
+	title = {{SoilGrids250m}: {Global} gridded soil information based on machine learning},
+	volume = {12},
+	issn = {1932-6203},
+	shorttitle = {{SoilGrids250m}},
+	url = {https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0169748},
+	doi = {10.1371/journal.pone.0169748},
+	abstract = {This paper describes the technical development and accuracy assessment of the most recent and improved version of the SoilGrids system at 250m resolution (June 2016 update). SoilGrids provides global predictions for standard numeric soil properties (organic carbon, bulk density, Cation Exchange Capacity (CEC), pH, soil texture fractions and coarse fragments) at seven standard depths (0, 5, 15, 30, 60, 100 and 200 cm), in addition to predictions of depth to bedrock and distribution of soil classes based on the World Reference Base (WRB) and USDA classification systems (ca. 280 raster layers in total). Predictions were based on ca. 150,000 soil profiles used for training and a stack of 158 remote sensing-based soil covariates (primarily derived from MODIS land products, SRTM DEM derivatives, climatic images and global landform and lithology maps), which were used to fit an ensemble of machine learning methods—random forest and gradient boosting and/or multinomial logistic regression—as implemented in the R packages ranger, xgboost, nnet and caret. The results of 10–fold cross-validation show that the ensemble models explain between 56\% (coarse fragments) and 83\% (pH) of variation with an overall average of 61\%. Improvements in the relative accuracy considering the amount of variation explained, in comparison to the previous version of SoilGrids at 1 km spatial resolution, range from 60 to 230\%. Improvements can be attributed to: (1) the use of machine learning instead of linear regression, (2) to considerable investments in preparing finer resolution covariate layers and (3) to insertion of additional soil profiles. Further development of SoilGrids could include refinement of methods to incorporate input uncertainties and derivation of posterior probability distributions (per pixel), and further automation of spatial modeling so that soil maps can be generated for potentially hundreds of soil variables. Another area of future research is the development of methods for multiscale merging of SoilGrids predictions with local and/or national gridded soil products (e.g. up to 50 m spatial resolution) so that increasingly more accurate, complete and consistent global soil information can be produced. SoilGrids are available under the Open Data Base License.},
+	language = {en},
+	number = {2},
+	urldate = {2023-07-05},
+	journal = {PLOS ONE},
+	author = {Hengl, Tomislav and Jesus, Jorge Mendes de and Heuvelink, Gerard B. M. and Gonzalez, Maria Ruiperez and Kilibarda, Milan and Blagotić, Aleksandar and Shangguan, Wei and Wright, Marvin N. and Geng, Xiaoyuan and Bauer-Marschallinger, Bernhard and Guevara, Mario Antonio and Vargas, Rodrigo and MacMillan, Robert A. and Batjes, Niels H. and Leenaars, Johan G. B. and Ribeiro, Eloi and Wheeler, Ichsani and Mantel, Stephan and Kempen, Bas},
+	month = feb,
+	year = {2017},
+	note = {Publisher: Public Library of Science},
+	keywords = {Agricultural soil science, Forecasting, Glaciers, Machine learning, Remote sensing, Shannon index, Soil pH, Trees},
+	pages = {e0169748},
+	file = {Full Text PDF:/Users/dorme/Zotero/storage/VS5FRVSK/Hengl et al. - 2017 - SoilGrids250m Global gridded soil information bas.pdf:application/pdf},
+}
+
+@article{davis:2017a,
+	title = {Simple process-led algorithms for simulating habitats ({SPLASH} v.1.0): robust indices of radiation, evapotranspiration and plant-available moisture},
+	volume = {10},
+	issn = {1991-959X},
+	shorttitle = {Simple process-led algorithms for simulating habitats ({SPLASH} v.1.0)},
+	url = {https://gmd.copernicus.org/articles/10/689/2017/},
+	doi = {10.5194/gmd-10-689-2017},
+	abstract = {Bioclimatic indices for use in studies of ecosystem function, species distribution, and vegetation dynamics under changing climate scenarios depend on estimates of surface fluxes and other quantities, such as radiation, evapotranspiration and soil moisture, for which direct observations are sparse. These quantities can be derived indirectly from meteorological variables, such as near-surface air temperature, precipitation and cloudiness. Here we present a consolidated set of simple process-led algorithms for simulating habitats (SPLASH) allowing robust approximations of key quantities at ecologically relevant timescales. We specify equations, derivations, simplifications, and assumptions for the estimation of daily and monthly quantities of top-of-the-atmosphere solar radiation, net surface radiation, photosynthetic photon flux density, evapotranspiration (potential, equilibrium, and actual), condensation, soil moisture, and runoff, based on analysis of their relationship to fundamental climatic drivers. The climatic drivers include a minimum of three meteorological inputs: precipitation, air temperature, and fraction of bright sunshine hours. Indices, such as the moisture index, the climatic water deficit, and the Priestley–Taylor coefficient, are also defined. The SPLASH code is transcribed in C++, FORTRAN, Python, and R. A total of 1 year of results are presented at the local and global scales to exemplify the spatiotemporal patterns of daily and monthly model outputs along with comparisons to other model results.},
+	language = {English},
+	number = {2},
+	urldate = {2023-07-05},
+	journal = {Geoscientific Model Development},
+	author = {Davis, Tyler W. and Prentice, I. Colin and Stocker, Benjamin D. and Thomas, Rebecca T. and Whitley, Rhys J. and Wang, Han and Evans, Bradley J. and Gallego-Sala, Angela V. and Sykes, Martin T. and Cramer, Wolfgang},
+	month = feb,
+	year = {2017},
+	note = {Publisher: Copernicus GmbH},
+	pages = {689--708},
+	file = {Full Text PDF:/Users/dorme/Zotero/storage/RUZWLGHT/Davis et al. - 2017 - Simple process-led algorithms for simulating habit.pdf:application/pdf},
+}

docs/source/users/pmodel/pmodel_details/pmodel_overview.md

@@ -67,7 +67,7 @@ The main steps are:
 ## Variable graph
 
 The graph below shows these broad model areas in terms of model inputs (blue) and
-modelled outputs (red) used in the P-model. Optional inputs and internal variables are
+modelled outputs (red) used in the P Model. Optional inputs and internal variables are
 shown with a dashed edge.
 
 ![pmodel.svg](pmodel.svg)

docs/source/users/pmodel/pmodel_details/soil_moisture.md

@@ -12,7 +12,6 @@ kernelspec:
   name: pyrealm_python3
 ---
 
-
 # Soil moisture effects
 
 At present, there are four approaches for incorporating soil moisture effects on
@@ -27,7 +26,9 @@ photosynthesis:
 * The experimental `rootzonestress` argument to {class}`~pyrealm.pmodel.pmodel.PModel`.
 * The `lavergne20_c3` and `lavergne20_c4` methods for
   {class}`~pyrealm.pmodel.pmodel.CalcOptimalChi`, which use an empirical model of the
-  change in $\beta$ with soil moisture.
+  change in the ratio of the photosynthetic costs of carboxilation and transpiration.
+  Altering this cost ratio - inconveniently also called $\beta$ - for soil moisture
+  stress provides a more complete picture of plant responses than GPP penalty factors.
 
 The first three of these approaches are described here, but see [here](optimal_chi) for
 details of the last method.
@@ -44,6 +45,14 @@ The factor requires estimates of:
 * a measure of local mean aridity ($\bar{\alpha}$, `meanalpha`), as the average annual
   ratio of AET to PET.
 
+```{admonition} Soil moisture
+The parameters used in the calculation of this factor were estimated using the
+plant-available soil water expressed as a fraction of available water holding capacity.
+That capacity was calculated for the observed data on a site by site basis using the
+`SoilGrids` dataset {cite:p}`hengl:2017a`. Ideally, soil moisture calculated in the same
+way should be used with this approach.
+```
+
 The functions to calculate $\beta(\theta)$ are based on four parameters, derived from
 experimental data and set in {class}`~pyrealm.constants.pmodel_const.PModelConst`:
 
@@ -54,7 +63,7 @@ experimental data and set in {class}`~pyrealm.constants.pmodel_const.PModelConst
 * An intercept (a, `soilmstress_a`) for the aridity sensitivity parameter $q$.
 * A slope (b, `soilmstress_b`) for the aridity sensitivity parameter $q$.
 
-The aridity measure (($\bar{\alpha}$) is first used to set an aridity sensitivity
+The aridity measure ($\bar{\alpha}$) is first used to set an aridity sensitivity
 parameter ($q$), which sets the speed with which $\beta(\theta) \to 0$ as $m_s$
 decreases.
 
@@ -64,7 +73,7 @@ $$
 
 Then, relative soil moisture ($m_s$) is used to calculate the soil moisture factor:
 
-$$`
+$$
     \beta(\theta) = q ( m_s - \theta^\ast) ^ 2  + 1
 $$
 
@@ -74,7 +83,7 @@ varies with changing soil moisture for some different values of mean aridity. In
 the examples below, the default $\theta_0 = 0$ has been changed to $\theta_0 =
 0.1$ to make the lower bound more obvious.
 
-```{code-cell} ipython3
+```{code-cell}
 :tags: [hide-input]
 
 from matplotlib import pyplot as plt
@@ -127,6 +136,8 @@ ax2.set_ylabel(r"Empirical soil moisture factor, $\beta(\theta)$")
 plt.show()
 ```
 
++++ {"user_expressions": []}
+
 ### Application of the {func}`~pyrealm.pmodel.functions.calc_soilmstress_stocker` factor
 
 The factor can be applied to the P Model by using
@@ -136,7 +147,7 @@ by the resulting factor. The example below shows how the predicted light use
 efficiency from the P Model changes across an aridity gradient both with and without the
 soil moisture factor.
 
-```{code-cell} ipython3
+```{code-cell}
 # Calculate the P Model in a constant environment
 tc = np.array([20] * 101)
 sm_gradient = np.linspace(0, 1.0, 101)
@@ -146,15 +157,26 @@ model = pmodel.PModel(env)
 model.estimate_productivity(fapar=1, ppfd=1000)
 
 # Calculate the soil moisture stress factor across a soil moisture gradient
-sm_stress = pmodel.calc_soilmstress_stocker(soilm=sm_gradient, meanalpha=0.5)
-penalised_gpp = model.gpp * sm_stress
+# at differing aridities
+
+gpp_stressed = {}
+
+for mean_alpha in [0.9, 0.5, 0.3, 0.1, 0.0]:
+    # Calculate the stress for this aridity
+    sm_stress = pmodel.calc_soilmstress_stocker(
+        soilm=soilm, meanalpha=mean_alpha, const=const
+    )
+    # Apply the penalty factor
+    gpp_stressed[mean_alpha] = model.gpp * sm_stress
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 :tags: [hide-input]
 
 plt.plot(sm_gradient, model.gpp, label="No soil moisture penalty")
-plt.plot(sm_gradient, penalised_gpp, label="Soil moisture penalty applied")
+
+for ky, val in gpp_stressed.items():
+    plt.plot(soilm, val, label=r"$\bar{{\alpha}}$ = {}".format(ky))
 
 plt.xlabel(r"Relative soil moisture, $m_s$, -")
 plt.ylabel(r"GPP")
@@ -183,9 +205,18 @@ The factor requires estimates of:
 
 * relative soil moisture ($m_s$, `soilm`), as the fraction of field capacity, and
 * a climatological estimate of local aridity index, typically calculated as total PET
-  over total precipitation for an appropriate period.
+  over total precipitation for an appropriate period, typically at least 20 years.
+
+```{admonition} Soil moisture
 
-The calculation of $f(\theta)$ is based on two functions of the aridity index: both
+The parameters used in the calculation of this factor were estimated using the
+plant-available soil water expressed as the ratio of millimeters of soil moisture over
+the total soil capacity. The soil water estimated using CRU data in SPLASH v1 model
+{cite:p}`davis:2017a`, which enforces a constant soil capacity of 150mm. Again, ideally,
+soil moisture calculated in the same way should be used with this approach.
+```
+
+The calculation of $\beta(\theta)$ is based on two functions of the aridity index: both
 power laws, constrained to take a maximum of 1 (no soil moisture stress penalty). The
 first function describes the maximal attainable level ($y$) and the second function
 describes a threshold ($\psi$) at which that level is reached. The parameters of these
@@ -199,7 +230,7 @@ y &= \min( a  \textrm{AI} ^ {b}, 1)\\
 \end{align*}
 $$
 
-```{code-cell} ipython3
+```{code-cell}
 from pyrealm.constants import PModelConst
 
 const=PModelConst()
@@ -235,7 +266,7 @@ $$
     \end{cases}
 $$
 
-```{code-cell} ipython3
+```{code-cell}
 # Calculate the soil moisture stress factor across a soil moisture
 # gradient for different aridity index values
 beta = {}
@@ -261,11 +292,10 @@ under drier conditions.
 
 As with  {func}`~pyrealm.pmodel.functions.calc_soilmstress_stocker`, the factor is first
 calculated and then applied to the GPP calculated for a model
-(~pyrealm.pmodel.pmodel.PModel.gpp). In the example below, the result is obviously just
-$\beta(\theta)$ from above scaled to the constant GPP.
-
-```{code-cell} ipython3
+({attr}`~pyrealm.pmodel.pmodel.PModel.gpp`). In the example below, the result is
+obviously just $\beta(\theta)$ from above scaled to the constant GPP.
 
+```{code-cell}
 for ai in ai_vals:
 
     plt.plot(sm_gradient, model.gpp * beta[ai], label= f"AI = {ai}")
@@ -281,7 +311,7 @@ plt.show()
 ```{warning}
 This approach is **an experimental feature** - see the
 {class}`~pyrealm.pmodel.pmodel.PModel` documentation. Essentially, the values for
-`rootzonestress` apply a penalty factor directly to $\beta$ in the calculation of
-optimal $\chi$. This factor is currently calculated externally to the `pyrealm` package
-and is not documented here.
+`rootzonestress` apply a penalty factor directly to the unit cost ratio $\beta$ in the
+calculation of optimal $\chi$. This factor is currently calculated externally to the
+`pyrealm` package and is not documented here.
 ```

pyrealm/pmodel/functions.py

@@ -527,7 +527,6 @@ def calc_soilmstress_stocker(
         \[
             \beta =
                 \begin{cases}
-                    0
                     q(m_s - \theta^{*})^2 + 1,  & \theta_0 < m_s <= \theta^{*} \\
                     1, &  \theta^{*} < m_s,
                 \end{cases}