Bioclimatic traits in statistical properties of daily photosynthetically active radiation

Abstract. In this paper, we present a methodology to analyze the stochastic component of daily solar radiation at the earth’s surface in the photosynthetically active spectral band. Extraterrestrial solar spectral irradiance from the SOLID project and in situ observed PAR from the FLUXNET data set are used to calculate daily time series of the clearness and clear-sky indices for 28 FLUXNET sites around the world for 1978-2014. We found that the shape of the probability distribution functions of the clearness and clear-sky indices exhibit a spatial pattern related to the Köppen climate classification and the Holdridge 5 life zones. According to the Köppen classification, oceanic, humid continental, and Mediterranean climates show bimodal distributions; semi-arid, temperate, subtropical, and desert climates show unimodal distributions with high dispersion; and tropical climates show unimodal distributions with low dispersion. Moreover, according to the Holdridge triangle, sites with bimodal distributions are concentrated in moist and wet forest life zones located in boreal and cool temperate regions and sub-humid and humid provinces. Unimodal distributions with high dispersion are concentrated in the moist forest life zone in 10 subtropical and tropical regions and humid province; and unimodal distributions with low dispersion are concentrated in dry forest, very dry forest, and thorn woodland in tropical and subtropical regions between arid and subhumid humidity provinces.


Data
Our data set consists of daily observations of incoming photosynthetic photon flux density (PPFD) and rainfall from 28 sites around the world from the FLUXNET data set (Baldocchi et al., 2001;Olson et al., 2004) (Fig. 1). Sites have different periods of record spanning from 1996 to 2014 and elevations from sea level to 1550 m (Table 1). We also use the Solar Spectral Irradiance (SSI) at the top of the atmosphere from the "First European Comprehensive Solar Irradiance Data Exploitation 5 project" (SOLID) (Haberreiter et al., 2017;Schöll et al., 2016) as input data for an atmospheric radiation transfer model (section 3). The SOLID spectral time series has a daily time resolution from 1978/7/11 to 2014/12/31 (13204 days) and covers the wavelength range between 0.5 and 1991.5 nm. Data from SOLID are available at http://projects.pmodwrc.ch/solid.
The wavelength domain for PPFD in the FLUXNET data set is 400-700 nm (Olson et al., 2004) and has units of µmol m −2 s −1 .
Information of Köppen climate classification and Holdridge life zones is taken from Rubel and Kottek (2010) and NEP-WCMC, respectively.

Methods
Daily radiation amount at a site on the earth's surface is the result of integrating instantaneous irradiance over the day length.
Surface instantaneous irradiance estimation comprises solar irradiance at the top of the atmosphere (TOA) and the physical properties of the atmosphere for the site and time of interest. We use SOLID data for TOA irradiance and the Beer-Lambert law to calculate light attenuation by the atmosphere. However, some atmospheric components as clouds, water vapor, and aerosols 5 are highly variable in space and time, which is troublesome when using the Beer-Lambert law. Therefore, we follow two approaches: 1) use the clear-day index (c) (also known as relative clearness index, clear day index, and normalized clearness index) to assess the effect of the variable components on total daily radiation, and 2) use the clearness index (k) to assess the whole atmospheric effect on total daily radiation. Both indices are defined in Eqs. (1), where P AR obs is the observed daily global PAR on a horizontal surface at the ground level, P AR 0 is the extraterrestrial daily global PAR on a horizontal surface, 10 and P AR cda is the daily global PAR on a horizontal surface on the ground for a cloudless, clean, and dry atmosphere. c = P AR obs P AR cda , k = P AR obs P AR 0 (1)

PAR at the top of the atmosphere
Photosynthetically active radiation at the top of the atmosphere (P AR 0 ) is calculated by integrating the daily SSI over the PAR spectral band (400-700 nm) to obtain the Total Spectral Irradiance (T SI P AR ), as shown in Eq.
transformations accounting for solar declination (δ) and latitude (φ), as shown in Eq.
where E 0 is the eccentricity correction factor of the earth's orbit and ω sr is the sunrise hour angle for the day.
3.2 Surface PAR for a cloudless, clean dry atmosphere 5 Surface PAR for a cloudless, clean, and dry atmosphere (P AR cda ) is the sum of the direct (P AR b ) and diffuse (P AR d ) components (P AR cda = P AR b + P AR d ). To calculate daily P AR b and P AR d on a horizontal surface at the ground level, we model the direct and diffuse instantaneous spectral irradiances and integrate them along the day length and the PAR spectral domain.
Following Iqbal (1983), we assume the cloudless, clean and, dry atmosphere to be composed by uniformly mixed gases (m) 10 and ozone (o). Using the Beer-Lambert law and integrating, daily P AR b is calculated as in Eq. (4).
where SSI 0,n,λ is the extraterrestrial spectral irradiance normal to the rays from the sun (obtained from SOLID), γ is the solar altitude varying from sunrise (sr) to sunset (ss), and τ ma,λ is the transmittance due to molecular absorbers of the atmosphere. 15 For the assumed atmosphere composition τ ma,λ = τ o · τ g , where τ o and τ g are the ozone and the mixed gases transmittance, respectively (see details in Iqbal, 1983, Sec.6.14). Assuming forward and backward scatterances of 0.5 and considering only the first pass of radiation through the atmosphere, P AR d can be calculated by where τ r,λ is the transmittance due to Rayleigh molecular scattering (see details in Iqbal, 1983, Sec.6.14).

20
Several atmospheric parameters are required by Eqs. (4) and ( where m r is relative air mass at standard pressure, m a is relative air mass at actual pressure, k o and k g are the absorption 5 attenuation coefficient for oxygen and mixed gases, and l o is the amount of ozone in cm (at normal temperature and pressure, NTP). We calculate k o,λ for any λ value using the Leckner (1978) interpolation of the classic Vigroux (1953) data. For calculating l o we interpolate Table 5.3.2 from Iqbal (1983), which is a reproduction of Robinson (1966, p.114). k o,λ is calculated by interpolating Table 6.13.1, which is a reproduction of Table 4 in Leckner (1978, p.146).

Statistical properties of k and c 10
We estimate the daily time series, annual cycles, autocorrelograms, and empirical probability density functions (PDF) of c and k. We separate the data of c and k by rainy and dry days, using precipitation as a proxy of cloudiness and water vapor in the atmosphere. Then, we inspect the seasonality of c and k by comparing the cumulative distribution function (CDF) of each month with the CDFs of the other months. The comparison of the CDFs is carried out visually and tested by using Kolmogorov-Smirnov (KS) and Anderson-Darling (AD) Goodness of Fit Tests (Dodge, 2008;Pearson, 1900;Scholz and 15 Stephens, 1987).

Results and Discussion
The time series, annual cycle, and autocorrelogram of PAR, c and k were calculated and plotted for each site. The results for AT-NEU are shown in Fig. 2, and those for all other sites are in Figs. S1 to S28. Notice that some values of c and k are greater than 1 during winter in sites with seasonal snow. This anomaly can be explained by the multiple reflection of light enhanced 20 by the snow cover.
The autocorrelograms of PAR, c, and k indicate a more marked annual cycle for PAR than for c and k (see The ACF of k is stronger than that of c in most sites ( Fig. 2, panels f and i). The calculation of k does not include any atmosphere, while c considers a cloudiness-sky, clean, and dry atmosphere. The high ACFs of k suggest that the atmosphere (specifically the air mass) has seasonality, that k does not manage to remove (Ianetz and Kudish, 2008).
ACFs show a period of 180 days approximately for all the sites studied, except in BR-Sa3 where it is 120 days (for both PAR and the indices), and in AT-Neu where it is 80 days (only for the indices). As seasonality depends largely on the movement 5 of the sun, the periods in sites nearest to the geographical Equator should be shorter than those of extratropics, as in BR-Sa3.
However, seasonality is also a function of local factors not dealt with in this paper. The differences in the periods of PAR and its indices in AT-Neu can be explained because the climatic and astronomical seasonalities are out of phase.
The pdfs of c and k reveal a certain degree of bimodality, or at least, some asymmetry in respect to the mean, as mentioned before by Hollands and Suehrcke (2013) is observed for k and c. This led us to suspect that k and c follow bimodal distributions in high latitudes, unimodal distributions with high dispersion in mid-latitudes, and unimodal distributions with low latitudes.
Looking for spatial patterns in the PDF shapes of daily c and k, we arranged sites following the Köppen classification and Holdridge life zones (Holdridge, 1947(Holdridge, , 1967. According to the Köppen classification, bimodality occurs in oceanic, humid continental, and Mediterranean climates; unimodal PDFs with high dispersion in tropical monsoon, tropical savanna, and

CH-Oe2
CH-Oe1  (1), the slope of lines (m) is given by P AR 0 /P AR cda . As P AR 0 does not consider the absorption by atmosphere components, it has a higher value than P AR cda that does it, being m always greater than 1. Sites in  Due to the multiple refractions of light by snow, we found values of c and k greater than 1 in sites where there is seasonal snow, during the periods in which it occurs.

US-NC2
The analysis of rainy and dry days revealed that c and k have similar statistical shape PDFs as indicated by Escobedo et al. (2009), but k exposed higher values than c. The differences between c and k are more accentuated in the extratropical northern hemisphere since the path length that energy must pass through varies considerably during the year. Besides, in the 5 extratropical northern hemisphere, c shows lower autocorrelations than k for all the lags analyzed, indicating that the air mass of the atmosphere has seasonality, and k does not manage to remove it.
By locating the 28 sites analyzed on the Holdridge life zones scheme (see Fig. 3), we noticed that the geographical location and climate greatly influence the statistical behavior of c and k. High latitudes sites exhibit bimodal distributions, arid to subhumid climates exhibit unimodal distributions with high dispersion, and humid tropical regions exhibit unimodal distributions 10 with low dispersion.
As rainfall is a proxy of cloudiness and water vapor, bimodal distributions vanish when data are divided into rainy and dry days, corroborating that each mode is related to clear and overcast skies conditions (Tovar-Pescador, 2008;Olseth and Skartveit, 1984).
On the seasonality, we found a clear definition of indices seasons at extratropical sites, unlike tropical sites where local 15 conditions play a determining role. Variations in the beginning and end of indices seasons do not show a notorious spatial