Model description

The Lagrangian stochastic SMOP model has been under development in the National Institute of Agronomic Research, France since 2003 (Jarosz et al., Reference Jarosz, Loubet and Huber2004). The model is designed to predict pollen dispersion and deposition downwind and within a flowering field. Displacement of individual pollen grains is given by the following two-dimensional stochastic differential equations.

(1)$${\rm d}u = a_u{\rm d}t + b_ud\xi _u{\rm d}x = u\,{\rm d}t$$

(2)$${\rm d}w = a_wd + b_wd\xi _w{\rm d}z = \lpar w-W_S\rpar \,{\rm d}t$$

where t is time; u and w are the horizontal and vertical air particle velocity components, respectively; dξ _u and dξ _w are random numbers drawn from Gaussian distributions with mean zero and variance dt; W_s is the gravitational settling velocity of pollen grains (terminal velocity of pollen grains in still air); and a_u, b_u, a_w and b_w are the Langevin coefficients (depending on horizontal and vertical averages of the air velocity components, U and W, respectively; the horizontal and vertical Eulerian velocity variances, $\sigma _u^2$ and $\sigma _w^2$, respectively; the shear stress, $\overline {u^{\prime}w^{\prime}}$; and the Lagrangian velocity time scale, T_L). All of these parameters are described in Loubet (Reference Loubet2000) and Loubet et al. (Reference Loubet, Cellier, Milford and Sutton2006).

The SMOP model computes pollen trajectories taking into account both canopy and pollen characteristics, and micrometeorological parameters (Loubet, Reference Loubet2000; Jarosz et al., Reference Jarosz, Loubet and Huber2004; Loubet et al., Reference Loubet, Cellier, Milford and Sutton2006). The mean wind speed profile which is implemented in the model above and in the canopy is defined, respectively, as (Raupach et al., Reference Raupach, Finnigan and Brunet1996) and (Garratt, Reference Garratt1992).

(3)$$u\,\lpar z\rpar = \left\{{\matrix{ {u\lpar h_c\rpar \times \exp \lpar \gamma \times {\textstyle{z \over {h_c}}}-1\rpar } \hfill & {z \lt h_r} \hfill \cr {{\textstyle{{u_\ast } \over k}}\,\left[{\ln \lpar {\textstyle{{z\,-\,d} \over {z_0}}}\rpar -\Psi_m\lpar {\textstyle{{z\,-\,d} \over L}}\rpar + \Psi_m\lpar {\textstyle{{z_0} \over L}}\rpar } \right]} \hfill & {z \ge h_r} \hfill \cr } } \right.$$

where h_c is the canopy height; u_* is the friction velocity; k is the von Karman constant; L is the Monin-Obukhov length; h_r is the height where the exponential and logarithmic profiles of the model join which is the height of the inflexion point of the wind speed profile; Ψ_m is the stability correction function given by Paulson (Reference Paulson1970) and Dyer (Reference Dyer1974); z ₀ = 0.1 × h_c is the roughness length and d = 0.7 × h_c is the displacement height (Kaimal and Finnigan, Reference Kaimal and Finnigan1994); and γ = 2 is the attenuation factor (Raupach et al., Reference Raupach, Finnigan and Brunet1996; Loubet et al., Reference Loubet, Cellier, Milford and Sutton2006). Required meteorological records were provided by an Agroclim station located in the experimental domain at Grignon, France. Data should be retrieved from the web site http://www6.paca.inra.fr/agroclim where global missions and objectives of Agroclim stations are given by the web site https://www6.paca.inrae.fr/agroclim/content/download/3214/32033/version/1/file/Plaquette+AGROCLIM+2020+FR_V1.pdf. For more information about the model, please visit the following website: https://www6.versailles-grignon.inra.fr/ecosys_eng/Staff/Photos-and-CV/L/Loubet-B/MODDAAS-SMOP.

Sensitivity analysis methods

Modelling a real-life process with a computer program, one is often faced with the problem of what values should be used for the input variables. Therefore, the output results have uncertainty associated with the input values. Stochastic approaches are usually used to obtain probabilistic output based on ‘randomly’ distributed parameters. Analysts desiring to conduct a sensitivity analysis using such an approach are obliged to select a set of values for their model inputs. This selection may be done in different ways: the random sampling is usually used in cases where the only general behaviour of flows are accounted for; the stratified sampling and the Latin hypercube sampling techniques are preferably chosen when studies need more accuracy: (i) what is the uncertainty in the output given the uncertainty in the input? and (ii) how important are the individual elements of the input with respect to the uncertainty in the output? (McKay et al., Reference McKay, Beckmann and Conover1979; Iman, Reference Iman1992). Latin hypercube sampling is a particular stratified sampling technique and in the absence of correlations between input variables, it generates a set of values by dividing the sample space of each variable into m intervals of equal probability (where m is the number of required model runs) and one value is randomly selected from each interval.

Regression analysis sensitivity

Using the least squared method, the model input variables X_j (1 ≤ j ≤ p) and the model output Y are used to construct a linear model which can be written as:

(4)$$Y_k = b_0 + \sum\limits_{\,j = 1}^p {b_j.X_{\,jk} + \varepsilon _k = {\hat{Y}}_k + \varepsilon _kk = 1\comma \;\ldots \comma \;m\comma \;} $$

where p is the number of the input variables; b_j (1 ≤ j ≤ p) is the ordinary regression coefficients (ORC) of the model output; ɛ_k is the regression residual; and $\hat{Y}_k$ is the regression-estimate of the model output for k^th model run which represents the best linear prediction. The signs of the coefficients b ₁, b ₂, …, b_p indicate whether Y_k increases (i.e., a positive coefficient) or decreases (i.e., a negative coefficient) as the corresponding X_jk values increases. The main drawback of the ORC method is the dependence of coefficients on physical units of variables and therefore these quantities cannot be compared with each other (Helton et al., Reference Helton, Johnson, Sallaberry and Storlie2006). Because of this, the ORC in Eqn (4) is expressed in the standardized regression coefficients (SRC) or normalized regression coefficients (NRC) forms (Saltelli et al., Reference Saltelli, Chan and Scott2000a) as:

$$\displaystyle{{\hat{Y}-\bar{Y}} \over {\hat{S}}} = \sum\limits_{\,j = 1}^p {b_j^{\lpar s\rpar } \cdot \displaystyle{{X_j-{\bar{X}}_j} \over {{\hat{S}}_j}}} \comma \;$$

and

$$\displaystyle{Y \over {\bar{Y}}} = \sum\limits_{\,j = 1}^p {b_j^{\lpar n\rpar } \cdot \displaystyle{{X_j} \over {{\bar{X}}_j}}} \comma \;$$

where $\lsqb \bar{Y}\comma \;\,\,\hat{S}\rsqb$ and $\lsqb \bar{X}_j\,\comma \;\,\hat{S}_j\rsqb$ are the means and standard deviations of (Y ₁,…, Y_m) and (X_{j 1},…, X_jm); and $b_j^{\lpar s\rpar }$ and $b_j^{\lpar n\rpar }$ are the SRC and NRC coefficients, respectively. These new coefficients may be obtained easily from the ORC:

$$b_j^{\lpar s\rpar } = b_j.\displaystyle{{{\hat{S}}_j} \over {\hat{S}}}\;{\rm and}\;b_j^{\lpar n\rpar } = b_j.\displaystyle{{{\overline X }_j} \over {\overline Y }}.$$

$b_j^{\lpar s\rpar }$ and $b_j^{\lpar n\rpar }$ measure respectively the effect of moving each input variable away from its mean by a fixed fraction of its standard deviation and the relative change $\Delta Y/\hat{Y}$ of Y due to relative change of $\Delta X_j/\hat{X}_j$, while the other variables remain constant. The regression study is usually associated with the coefficient of determination R ² which is computed using the variance of the predicted values in relation to the variance of the observed ones. More specifically, if R ² is close to 1 then the linear relationship is strong. The usefulness of the coefficient R ² in sensitivity analysis is limited by the fact that when additional variables are added to the regression equation, the value of R ² increases even when these variables do not significantly improve the regression equation (Janssen, Reference Janssen1994). Because of this, the adjusted $R_{adj}^2$ is usually used instead:

(5)$$R_{adj}^2 = 1-\lpar 1-R^2\rpar \times \displaystyle{{m-1} \over {m-p-1}}\comma \;$$

where m is the model runs and p is the number of the input variables. The adjusted $R_{adj}^2$ provides a measure of the extent to which the regression model can match the observed data and it is more appropriate for a multiple linear regression than the coefficient of determination R ² (Janssen et al., Reference Janssen, Heuberger and Sanders1992; Janssen, Reference Janssen1994).

Linear correlation analysis sensitivity

The linear correlation coefficient (LCC) measures the strength of the linear relationship between the input variable X_j and the output Y; it can be written as:

(6)$$\rho \lpar X_j\comma \;Y\rpar = \displaystyle{{\sum\nolimits_{k = 1}^m {\lpar X_{\,jk}-{\bar{X}}_j\rpar \lpar Y_k-\bar{Y}\rpar } } \over {\sqrt {\sum\nolimits_{k = 1}^m {{\lpar X_{\,jk}-{\bar{X}}_j\rpar }^2\sum\nolimits_{k = 1}^m {{\lpar Y_k-\bar{Y}\rpar }^2} } } }}\comma \;$$

where $\bar{X}_j = \sum\nolimits_{k = 1}^m {X_{jk}/m}$ and $\bar{Y} = \sum\nolimits_{k = 1}^m {Y_k/m}$. LCC is always ranged between −1 and +1, values close to either these extremes indicate a strong correlation, while small absolute values indicate little or no correlation.

Statistical performance indices

The quantitative evaluation of models is usually made through statistical performance analysis that is presented and widely used in many works. The statistical indices which we will use here are defined as (Hanna et al., Reference Hanna, Chang and Strimaitis1993; Chang & Hanna, Reference Chang and Hanna2004):

the fractional bias (FB) where the general expression is given by

(7)$${\rm FB} = 2 \times \displaystyle{{{\overline C }_p-{\overline C }_o} \over {{\overline C }_p + {\overline C }_o}}\comma \;$$

where C_p and C_o are model predictions and observations of pollen concentration; $\overline C$ is the average over the dataset. FB produces values ranged along a linear scale and the systematic bias refers to the arithmetic difference between C_p and C_o, advantages of FB is that it is a dimensionless number, symmetrical and bounded; the fractional bias values are ranged between −2.0 (extreme under-prediction) to +2.0 (extreme over-prediction).

Like FB, the geometric mean bias (MG) is a measure of mean bias and indicate only systematic errors which lead to always overestimation or underestimation of the measured values. MG is defined as

(8)$${\rm MG} = \exp \lpar {\overline {\ln \lpar C_p\rpar } - {\overline {\ln \lpar C_o\rpar } } \rpar \comma \;} $$

It is clear to see that MG value greater than 1 implies that the model underestimates and an MG value less than 1 that the model overestimates.

The normalized mean square error (NMSE) and the geometric mean variance (VG) are measures of scatter and reflect both systematic and unsystematic (random) errors. NMSE is an estimator of the overall deviations between predicted and measured values.

(9)$${\rm NMSE} = \displaystyle{{\overline {{\lpar {C_o-C_p} \rpar }^2} } \over {\overline {C_o} \times \overline {C_p} }}\comma \;$$

The normalization by the product $\overline C _o \times \overline C _p$ ensures that the NMSE will not be biased towards models that over-predict or under-predict. Smaller values of NMSE denote better model performance, however high NMSE values do not necessarily mean that a model is completely wrong. The VG is given by

(10)$${\rm VG} = \exp \overline {\lpar {{\lpar {\ln \lpar C_o\rpar -\ln \lpar C_p\rpar } \rpar }^2} \rpar } $$

According to Ahuja and Kumar (Reference Ahuja and Kumar1996) and Kumar et al. (Reference Kumar, Luo and Bennett1993), the performance of a model can be deemed acceptable if: |FB| < 0.30, 0.75 < MG < 1.25, 1.00 ≤ VG < 4.00 and NMSE < 1.50.

A perfect model would have MG and VG = 1; and FB, and NMSE = 0. However, MG and VG are known to be strongly influenced by extremely high and low values, while FB and NMSE are more influenced by high observed and predicted concentrations (Chang and Hanna, Reference Chang and Hanna2004).

Methodology: sensitivity analysis and inferred pollen emission rates

The main objective of this paper is to infer the pollen emission rate from measurements and simulations of pollen deposition rates in a canopy gap for hybrid maize of 2.02 m height. The experiment was carried out unfortunately only over 1 day on 8 August 2008 on a 0.5 ha (120.4 × 41 m²) maize crop in Grignon, France (48°50′N, 1°56′E, 101 m a.s.l.) with a row spacing of 0.8 m and plant density of 9.5 plants/m², however, the developed methodology should be applied for multiple days. Measurements were done when released pollens by plants were high and the leaf area density is expressed in horizontal and vertical projections where ranges of variability are defined below. The ground deposition rate was recorded in a canopy gap of 2.4 m diameter at seven distances (0.0, 0.4, 0.8, 1.2, 1.6, 2.0 and 2.4 m) on a line parallel to the row and to the wind direction, starting where the airflow entered the gap and ending at the downwind edge of the plants encircling the gap. A container was placed at the height of the tassel base in the first inter row of maize upwind of the canopy gap to measure pollen emission rate. These measurements will be considered to provide good approximations of the real pollen emission rates. Measurements were recorded during 30 min intervals with six intervals recorded between 10:15 and 13:15 UTC where the maximum pollen emission rate typically occurs (Scott, Reference Scott1970; Gregory, Reference Gregory1973).

Deposition rates were measured using containers (diameter = 50 mm, height = 70 mm) filled with 30 ml of an electrolyte solution (Coulter Isoton, Beckman USA) and the amount of pollen in the containers were counted by using an automatic counter (Coulter Multisizer III, Beckman, USA). The wind speed and global radiation were recorded at 5 m above the canopy while relative humidity, rainfall, water vapour pressure, water vapour pressure deficit, presence of wetness and temperature were recorded at plant height. Each meteorological variable was averaged over successive 30 min intervals, as for the deposition measurements.

To start, sensitive analysis of the pollen deposition in the canopy gap according to seven input variables of SMOP was performed. These variables describe pollen and canopy characteristics as well as pollen interactions with its environment. There are four continuous variables: two variables of the Gaussian distribution of pollen settling velocity (mean, meanW_s, ranged from 0.15 m/s, which corresponds to dry pollen to 0.35 m/s, which corresponds to fresh pollen and standard deviation, StdW_S, ranged from 0.01 to 0.06 m/s) (Aylor, Reference Aylor2002; Loubet et al., Reference Loubet, Jarosz, Saint-Jean and Huber2007; Marceau et al., Reference Marceau, Saint-Jean, Loubet and Huber2012); horizontal and vertical projections of the leaf area density LAD_x and LAD_z, respectively (LAD_x ranged from 0 to 0.55 m²/m³ and LAD_z ranged from 0 to 0.65 m²/m³, Drouet (Reference Drouet2003)). More specifically, LAD_x influences deposition due to the gravitational settling, and for example, LAD_x = 0 m²/m³ means that no deposition occurs within the canopy due to the gravitational settling. While LAD_z influences deposition due to the inertial impaction, LAD_z = 0 m²/m³ means that no deposition occurs within the canopy due to the inertial impaction. In addition, there are three 2-level variables (Yes or No): deposition on leaves (D_f); resuspension of pollen on vegetation and ground (R_p); and horizontal pollen fluctuations due to the wind (1D). In case when no deposition on leaves will be considered (D_f = 1) we impose the probability of deposition trajectories equal to zero. A detailed description of the three 2-level variables D_f, R_p and 1D as well as a description of their eight possible combinations are given in Table 1.

Table 1. Description of the 2-level parameters which constitute the eight interactions used to conduct a sensitivity analysis of the model with a description of each of the eight possible combinations of the three variables D_f, R_p and 1D

Due to the difficulty to describe correctly the turbulence flow in the canopy gap, as we detail in the Discussion section, the following SMOP parameters were estimated using meteorological data; wind speed (4.06 m/s), canopy height (2.02 m), $\sigma _u/u_\ast$ and $\sigma _w/u_\ast$ above the canopy were set to 2.6 and 1.25, respectively, neutral conditions were assumed. The crop length upwind and downwind of the gap (105 and 13 m) and the assumed released pollen rate to the atmosphere was set at 100 grains/m²/s. Figure 1 describes the geometric configuration of the experimental study.

Fig. 1. Configuration of the computational domain with the mean wind direction and canopy dimensions. Seven containers were placed in the canopy gap at ground level, on a line parallel to the wind direction. Another container was placed in the first inter-row of maize, upwind of the canopy gap and at the height of the tassel base.

The output of the model considered here is the pollen deposition rate in the centre of the canopy gap which is a typical practice configuration for measurements of pollen deposition on the ground required to infer the pollen emission rate. To do this and to examine how the uncertainty in the variable estimates affects the deposition rate, the Latin hypercube sampling is used. This technique, as noted above, allows exploration of the full variable space with a greatly reduced number of model runs (Iman, Reference Iman1992). Therefore, only 105 sets of variables are created to run the model.

In the first approach, the sensitivity methods described in the previous subsections were applied independently to the seven variables. As a second approach, the three 2-level variables were replaced by a unique discrete variable that has been named hereafter as interaction representing the eight combinations of the three 2-level variables. To finish, for a set of SMOP parameters, inferences of pollen emission rate were obtained corresponding to the time intervals where pollen deposition rates were measured. Pollen emission rate was inferred by assuming that the ratio between the simulated pollen deposition rate (D_simulated) and the measured pollen deposition rate (D_real) is equal to the ratio between the simulated pollen release rate (S_simulated = 100 grains/m²/s) and the inferred pollen emission rate (S_real).

(11)$$\displaystyle{{S_{simulated}} \over {S_{real}}} = \displaystyle{{D_{simulated}} \over {D_{real}}}\Leftrightarrow S_{real} = S_{simulated} \times \displaystyle{{D_{real}} \over {D_{simulated}}}$$

For a given time interval, the pollen release rate is inferred from the measured deposition rate at seven distances in the canopy gap following Eqn (11). These inferred estimations are used to compute an average and a standard error of pollen emission rate. To assess the influence of SMOP input variables on the inferred pollen release rate, SA is done using Eqn (11) for different interactions.