Direct effect of aerosols on solar radiation and gross primary production in boreal and hemiboreal forests

2.1 Sites and data sets

We used data from five remote forest sites located at the middle and relatively high latitudes in Finland, Estonia and Russia: SMEAR I (Värriö, Finland), SMEAR II (Hyytiälä, Finland), SMEAR Estonia (Järvselja, Estonia), Zotino (Zotino, Krasnoyarsk region, Russia) and Fonovaya (Tomsk region, Russia). Figure 1 illustrates the sites' locations on the map, while Table 1 summarizes information on the data sets used in the present study. The information regarding the instruments used can be found in separate papers that describe the stations (listed below) and on the AVAA Internet site for SMEAR I and II (https://avaa.tdata.fi/web/smart, last access: 30 June 2018).

Figure 1Location of the sites (see Table 1 for latitudes and longitudes of the different stations).

Download

SMEAR II (Hari and Kulmala, 2005) is located at Hyytiälä Forestry Field Station in a conifer boreal forest near a lake in central Finland. The site is a managed, 55-year old Scots pine (Pinus sylvestris L.) forest stand with a closed canopy and an average tree height of ca. 18 m.

SMEAR I (Hari et al., 1994) is the site characterized by the highest latitude and is located in Finnish Lapland. The site is comprised of ca. 60-year old Scots pines with a rather open canopy. The average tree height is ca. 10 m.

SMEAR Estonia (Noe et al., 2015) is located in a hemiboreal forest zone and the stands in the tower footprint consist of a mixture of coniferous, Scots pine and Norway spruce (Picea abies L. Karst), and deciduous, silver birch (Betula pendula Roth.) and downy birch (Betula pubescens Ehrh.), trees. Because of the high heterogeneity and the mosaic of stands of hemiboreal forests the stand age is 65-years old on average ranging from 43-year old larch stands to 120-year old pine stands. The average ages of the dominating species are 102 years for pine, 79 years for spruce and 68 years for birch. Therefore, the tree height is variable: 22 m on average (ranging between 10 and 30 m).

Fonovaya in the Tomsk region is the most southern site (Matvienko et al., 2015). It is located on the Ob River in western Siberia, Russia, and forest stand is represented by mixed forest. The stand is dominated by 50-year old Scots pine, 45-year old birch (Betula verrucósa) and 27-year old aspen (Populus tremula). The average tree height is ca. 30 m, ranging from 25 m for birch up to 40 m for pines.

Zotino (Heimann et al., 2014; Park et al., 2018) is located near the Yenisei River in central Siberia. The forest is dominated by a 100-year old Scots pine (Pinus sylvestris L.) forest stand with an open canopy and an average canopy height of ca. 20 m.

Thus, we have data sets from five sites with three of them represented by pine stands and two represented by mixed forests; all of the sites are located at mid-latitudes. Note that currently the data from these five stations form the largest possible set of simultaneous atmospheric observations on trace gases, meteorology, solar radiation and aerosols, conducted in boreal and hemiboreal forests in Eurasia. Some of these sites lack certain components necessary for the analysis and therefore we use parameterizations when needed (and possible). For example, the diffuse fraction of solar radiation, fdif_bb, was parameterized in the data set from Fonovaya. Following Gu et al. (2002) and Alton (2008), we used the formula

where R_d is diffuse radiation, R_g is global radiation and $x=R_{\mathrm{g}}/R_{\mathrm{TOA}}$ is the ratio of the measured global radiation to the modelled radiation at the top of atmosphere. For x<0.28 we used fdif_bb=0.95, and for x>0.75 we used fdif_bb=0.1. The radiation at the top of atmosphere was calculated as R_TOA=I₀cos(sza), where I₀ is the solar constant and sza is the solar zenith angle.

The data used in this study correspond to the peak growing season defined as the time period with maximum GPP. Typically it includes June–August and part of May and September for all the sites except SMEAR I, where it includes July–August and part of June and September. For consistency, we used June–August data for SMEAR II, SMEAR Estonia, Fonovaya and Zotino and July–August data for SMEAR I. The analysis was performed for daytime (09:00–15:00 local time), 30 min averaged data.

The aerosol number-size distribution at Fonovaya was measured using two instruments, which do not overlap with respect to the size range. Therefore, there is a gap between 200 and 300 nm in the distribution. This gap was filled with an average value between the two adjacent points (one point was added between the two points). Apart from that, we did not apply gap-filling for aerosol data in this study, and used only good quality data. We did not fill gaps in the solar radiation data.

Note, that the CO₂ flux was obtained by the micrometeorological (eddy covariance) method at the SMEAR stations and Zotino, while the gradient method was applied to obtain the CO₂ flux from the Fonovaya data set. GPP was calculated using the following formula:

where TER is the total ecosystem respiration and NEE is the net ecosystem exchange. TER was obtained for all ecosystems using the nighttime method of CO₂ flux partitioning (Reichstein et al., 2005). The partitioning method at SMEAR I and II was based on soil temperature which makes the method more precise (Kolari et al., 2009; Lasslop et al., 2012), while at other sites it was based on air temperature (soil temperature is not measured). According to Lasslop et al. (2012), a possible consequence of an average GPP calculated during daytime (09:00–15:00), is a decrease in GPP calculated using soil temperature of about 0.5 µmol s⁻¹ m⁻² compared with GPP calculated using air temperature. In addition, in the data sets from the SMEAR stations absent GPP points were modelled (Kolari et al., 2009), while in the data sets from the Russian stations we did not fill the gaps. For the Zotino data set a good data percentage was 55 % on average (Park et al., 2018), and for the Fonovaya data set it was approximately 30 %.

Studying aerosol–radiation interactions, we used the Solis clear-sky model (see more in Sect. 2.2). The input parameters for this model are the aerosol optical depth at 700 nm (AOD₇₀₀) and precipitable water (PW). We used AOD at 675 nm and PW from AERONET (Holben et al., 1998); in particular, data from the following AERONET sites were used: Hyytiala (for SMEAR II), Sodankyla (for SMEAR I), Tomsk22 (for Fonovaya) and Toravere (for SMEAR Estonia). The AERONET sites are in the immediate vicinity of the Fonovaya and SMEAR II stations, SMEAR Estonia is 50 km from Toravere and SMEAR I is approximately 70 km from Sodankyla. We used Version 2, Level 2 data (cloud screened and quality controlled), except for Tomsk22 where we used Version 3, Level 2 data. The input data were averaged over daytime.

Currently aerosol characteristics are not measured at Zotino and there are no AERONET sites nearby; therefore aerosol–radiation interactions have not been studied for Zotino. However, radiation and CO₂ flux, which are necessary for the radiation–photosynthesis analysis, are measured; consequently, the Zotino data set has been included in this study.

2.2 The Solis clear-sky model

To distinguish between the effects of aerosols and clouds on solar radiation, we used a simplified broadband version of a clear-sky radiative transfer model (RTM) – Solis (Ineichen, 2008). This is a quasi-physical model which means that it employs Lambert–Beer relations for the general estimates of irradiance while using parameterizations for the total optical depths. The input parameters are the aerosol optical depth at 700 nm (AOD₇₀₀) and precipitable water (PW). Parameterizations are developed for the following range of parameters: sea level < altitude < 7000 m, 0 < AOD₇₀₀ < 0.45 and 0.2 cm < PW < 10 cm.

Solis parameterizations were obtained using the “urban” type of aerosol size distribution in the full RTM (Ineichen, 2008). The difference between calculations for “urban” and “rural” types of aerosol for the same AOD₇₀₀=0.18 was reported by Mueller et al. (2004). This AOD₇₀₀ value is larger than the typical values for all of the places that we consider in this study; hence, we expect smaller errors due to the inconsistent aerosol type. The difference reported by Mueller et al. (2004) was negligible for direct irradiance, whereas global irradiance for “rural” aerosol was around 20 W m⁻² larger during the daytime. Given that for clear-sky conditions between 09:00 and 15:00 and during the growing season, global irradiance drops to ∼600 W m⁻², this difference introduces a maximum error of 3 %. More tests for several sites in the USA and comparison between Solis and two more simplified clear-sky models, Bird and REST2, were reported by Sengupta and Gotseff (2013). Solis is the optimal model from the point of view of the required parameters: it performs only slightly worse than the other two models, and it requires only two input parameters.

We used Solis to model both global and diffuse irradiance. The horizontal global irradiance, I_gh, was calculated as follows:

while the direct irradiance, I_dir, was calculated as

Here $I_{\mathrm{0}}^{\prime }$ is a common modified (enhanced) irradiance defined in Ineichen (2008), τ_b and τ_g are the direct and global total optical depths, respectively, b and g are the fitting parameters, and sza is the solar zenith angle. Diffuse radiation can be found from the global radiation balance using the following equation:

All of the parameterizations used in the model are given in Ineichen (2008). The algorithm used for the calculations is written in Fortran. For the calculation of the sza, the online calculator Solar Position Algorithm (SPA) was used (available from http://www.nrel.gov, last access: 22 February 2018).

2.3 Accounting for the difference between PAR and broadband radiation

Note that the measurement methods of solar radiation are different at the different sites. At SMEAR II before 2009 and Fonovaya only broadband radiation is measured, while at SMEAR I and Zotino only photosynthetically active radiation (PAR) is measured. At SMEAR II after 2009 both global PAR and broadband global radiation, as well as diffuse PAR are measured. Broadband shortwave radiation, referred to as broadband radiation, is the radiation in the spectral range between 0.3 and 4.8 µm, while PAR is the radiation relevant for photosynthesis, i.e. in the range of wavelengths between 400 and 700 nm. The former is typically measured using thermopile pyranometers (energy sensors) and reported in energy units (W m⁻²), while the latter is measured using quantum sensors and is reported in quantum units (µmol s⁻¹ m⁻²).In the following, we quantify the ratio between global PAR and global broadband radiation, and the ratio between the diffuse fraction of PAR and the diffuse fraction of broadband radiation.

Following Ross and Sulev (2000), the ratio of PAR (µmol s⁻¹ m⁻²) to the broadband radiation (W m⁻²) is called the PAR quantum efficiency χ. By dividing global PAR by broadband global radiation at SMEAR II and finding its mean value (the values were obtained for the growing season and years listed in Table 1), χ_glob=2.06 µmol s⁻¹ W⁻¹ was obtained, and was somewhat higher than 1.81 µmol s⁻¹ W⁻¹ reported by Ross and Sulev (2000) for Estonia.

We explain the potential difference between the diffuse fraction of PAR and the diffuse fraction of broadband radiation as follows. Aerosol particles influence a certain part of the spectra of solar irradiance depending on the particle size distribution (Seinfeld and Pandis, 2016). This effect can be better understood using a dimensionless size parameter πd_p∕λ, where λ is the wavelength of incident irradiance and d_p is the particle diameter. If $\mathit{\pi }{d}_{\mathrm{p}}/\mathit{\lambda }\ll \mathrm{1}$, then Rayleigh scattering is prevailing, while $\mathit{\pi }{d}_{\mathrm{p}}/\mathit{\lambda }\gg \mathrm{1}$ means that the scattering properties of the particles are determined by the geometrical optics, i.e., so-called “geometric scattering”. The characteristics of the aerosol distribution become important for solar irradiance if $\mathit{\pi }{d}_{\mathrm{p}}/\mathit{\lambda }\sim \mathrm{1}$. For boreal forests where the particle size distribution is typically governed by the well-pronounced mode with the geometric mean diameter d_p≈100 nm (Dal Maso et al., 2008), the effective interaction of aerosols with solar radiation occurs in the ultraviolet range (100–400 nm) and in the blue part of the optical spectrum (400–500 nm). Compared to PAR (400–700 nm), the essential part of the broadband radiation energy is contained in the near infrared part of the spectrum (700–1400 nm), which is not affected by the relatively small particles prevailing in aerosol distributions typical for boreal forest. In other words, the effect of aerosols on the diffuse fraction of solar irradiance can be different for PAR and broadband radiation. Qualitatively, one would expect that both the increase in diffuse radiation and the decrease in global radiation would be more pronounced for PAR, meaning that the diffuse fraction of PAR would be higher under the same aerosol conditions. In order to make our analysis consistent and to be able to compare results from different sites, we performed an analysis of the data from SMEAR Estonia to compare the diffuse fractions of PAR and broadband radiation. The data set includes 4 years, from 2014 to 2017.

The measurements at SMEAR Estonia are made using an energy sensor. The hyperspectral radiometer SkySpec (Kuusk and Kuusk, 2018) is a purpose-built instrument for the automated measurement of global and diffuse sky irradiance. To obtain PAR radiation, the spectral data are converted from energy units to quantum units and are integrated over a 400–700 nm spectral range. Integration over the whole available spectral range in energy units is used for simulating a thermopile pyranometer.

2.4 The condensation sink as a measure of aerosol loading

As previously mentioned in Sect. 1, this study can be considered as a part of the project quantifying the COBACC feedback loop using ground-based measurements. The condensation sink (CS) is a typical parameter calculated from a ground-based aerosol number-size distribution and characterizing aerosol loading. It can be related to measured organic vapours, making it possible to study the effect of the formation and growth of secondary aerosol for photosynthesis. In addition, the CS was chosen in the original study of the COBACC feedback loop by Kulmala et al. (2014); therefore, for the sake of comparison we resort to this parameter.

The CS is calculated from the particle number-size distribution as follows: (Kulmala et al., 2001):

where D_v is the diffusion coefficient of the condensing vapour, n(d_p) is the particle number-size distribution, and β is the Fuchs–Sutugin coefficient. The physical meaning of the CS is the inverse time needed for vapours to be taken up by existing aerosol particles. Similar to the scattering coefficient and the AOD, the value of the CS depends on aerosol surface area. However, the contribution of large particles to the CS diminishes with increasing particles diameter, in contrast to the aerosol surface area. This effect becomes pronounced for particles with a diameter larger than about 300 nm. For boreal forests, where the mode with a geometric mean diameter of around 100 nm dominates the particle number-size distribution, the CS can be assumed to be directly proportional to the aerosol surface area (Ezhova et al., 2018). Thus, one can assume that the CS is an appropriate measure of atmospheric aerosols for the radiation studies in boreal forests. The connection between this parameter and AOD₅₀₀ for boreal forests is discussed in more detail in Appendix A.

2.5 The LUE and PAR analysis to assess the effect of diffuse radiation on GPP

There is a strong evidence that GPP dependence on R_d∕R_g is non-linear and has a maximum (Alton, 2008; Moffat et al., 2010; Park et al., 2018); however, this maximum is not always well pronounced. In what follows we explain the GPP maximum based on the ecosystem LUE and PAR dependences on R_d∕R_g. Following Cheng et al. (2016), we defined the LUE as the GPP per unit PAR, therefore $\mathrm{GPP}=\mathrm{LUE}\cdot \mathrm{PAR}$. Strictly speaking, the LUE is defined as the GPP per unit absorbed PAR, i.e. ${\mathrm{PAR}}_{\mathrm{abs}}=f\mathrm{APAR}\cdot \mathrm{PAR}$, where fAPAR is the fraction of the absorbed PAR. The fraction of the absorbed PAR depends on the leaf area index (LAI), the solar zenith angle and other factors. This dependence for boreal forests was studied in Maijasalmi (2015) and Hovi et al. (2016). Based on results reported by Hovi et al. (2016), fAPAR for tree heights larger than 10 m and at a moderate zenith angle (40–60^∘) can be estimated to be between 0.8 and 0.9. One can obtain the LUE defined by the absorbed PAR by dividing the LUE used in the present study by 0.8–0.9.

For all ecosystems with a sufficiently deep canopy and a high leaf area index, the LUE is expected to increase with R_d∕R_g, as a larger fraction of available photons can penetrate inside the canopy and can be used for photosynthesis. Some studies (Niyogi et al., 2004) have reported a linearly growing dependence of the LUE on R_d∕R_g. Furthermore, a decrease of PAR with R_d∕R_g can be expected, as an increase in the diffuse fraction of global irradiance corresponds to the enhancement of the scattering and reflecting properties of the atmosphere due to the presence of aerosols or clouds. Therefore, for each site the dependence of the LUE and PAR on R_d∕R_g was investigated separately, after which the GPP dependence on R_d∕R_g was derived from these two dependences.

Again, in order to have consistent data sets, we recalculated R_d∕R_g obtained from the PAR measurements so that it existed in terms of broadband radiation at all the sites. Conversely, broadband global radiation from Fonovaya and part of the SMEAR II data set was recalculated to PAR when investigating the LUE and PAR dependence on R_d∕R_g. We multiplied the global radiation by χ_glob=2.06 µmol s⁻¹ W⁻¹ in order to get PAR in quantum units. The PAR quantum efficiency was chosen to be equal to the one at SMEAR II, as for the daytime and similar solar zenith angles it is mostly aerosol dependent; SMEAR II and Fonovaya generally have similar aerosol loading values which can be confirmed by their similar CS values (Fig. 5). Considering the LUE, we filter out the data with low global irradiance, R_g⩽100 W m⁻² (PAR ⩽ 200 µmol s⁻¹ m⁻²). Below this critical R_g, the LUE shows significant scatter (it is high for the low radiation values); therefore we excluded these data from analysis.

3.1 Aerosol effect on solar radiation

3.1.1 Criterion of clear sky based on the Solis clear-sky model

To understand the importance of the clear-sky criterion for the diffuse fraction of global radiation, we report the model test against diffuse and global irradiance measurements at SMEAR Estonia ( the Solis clear-sky model is described in Sect. 2.2). Examples of the global and diffuse radiation diurnal cycles for clear and cloudy days are displayed in Fig. 2. Note that on a cloudy day with patchy clouds, the AOD can still be measured. Model results for a cloudy day report the global and diffuse clear-sky radiation for AOD and PW measured on that day. In general, the model performs well during clear-sky conditions as can be seen in Fig. 2a. This is in accordance with the results of Sengupta and Gotseff (2013) who reported good performance of the model for clear-sky conditions at several sites in the USA. On a cloudy day (Fig. 2b), there were times (e.g. at 09:00 and 17:00) when the measured and modelled global irradiance (R_g,meas and R_g,mod) were nearly equal while the measured diffuse irradiance (R_d,meas) was significantly higher than the modelled one (R_d,mod). The criterion of clear sky based on the comparison between the modelled and measured global irradiance, i.e. involving only global irradiance (Kulmala et al., 2014), does not filter out these points.

Figure 3The modelled vs. measured diffuse fraction of global irradiance (SMEAR Estonia, 2016). Closed symbols represent the criterion of clear sky based on diffuse and global irradiance ($R_{\mathrm{g},\mathrm{meas}}/R_{\mathrm{g},\mathrm{mod}}>\mathrm{0.9}$, $R_{\mathrm{d},\mathrm{meas}}/R_{\mathrm{d},\mathrm{mod}}>\mathrm{0.8}$), the closed symbols and the open symbols combined represent the criterion of clear sky based on global irradiance ($R_{\mathrm{g},\mathrm{meas}}/R_{\mathrm{g},\mathrm{mod}}>\mathrm{0.9}$). The solid line illustrates the ideal 1:1 ratio between the modelled and measured data sets.

Download

A further illustration of the simplified criterion and its consequences for the diffuse fraction of global irradiance is given in Fig. 3, displaying the diffuse fraction of global irradiance under clear-sky conditions at SMEAR Estonia (summer 2016). The open symbols in combination with the closed symbols correspond to the data obtained using the clear-sky criterion involving global radiation alone:

$$\begin{array}{}\text{(7)}& R_{\mathrm{g},\mathrm{meas}}/R_{\mathrm{g},\mathrm{mod}}\ge \mathrm{0.9}.\end{array}$$

The closed symbols, in comparison, show the data obtained with the criterion involving both the diffuse and global radiation:

$$\begin{array}{}\text{(8)}& R_{\mathrm{g},\mathrm{meas}}/R_{\mathrm{g},\mathrm{mod}}\ge \mathrm{0.9}\mathrm{and}\mathrm{0.8}\le R_{\mathrm{d},\mathrm{meas}}/R_{\mathrm{d},\mathrm{mod}}\le \mathrm{1.2}.\end{array}$$

The criterion of clear sky suggested here is based on the results of Sengupta and Gotseff (2013), who determined a rms error of the linear regression corresponding to the measured global radiation and modelled global radiation using Solis for eight sites. This rms error did not exceed 10 % (cf. Equation 7 using global radiation), and was generally lower than that. As for diffuse radiation, we chose a 20 % difference between measured and modelled diffuse radiation in Equation (8). This difference was chosen based on the estimated 18 W m⁻² error in diffuse radiation between the full Solis radiative transfer model and measurements, reported by Ineichen (2008), by assuming a typical diffuse radiation value of 120 W m⁻², and by adding a 5 % error between the simplified model and the full Solis radiative transfer model. The closed symbols in Fig. 3 show a good agreement between the measured and modelled diffuse fractions of solar irradiance R_d∕R_g, while a considerable portion of the open symbols has large measured R_d∕R_g – meaning that they contain a large fraction of cloud-influenced data. Therefore, the criterion of clear sky based on the global and diffuse radiation can be used to detect clear-sky data when it is important to separate the effect of aerosol and clouds on diffuse radiation.

Note that the diffuse fraction of solar irradiance due to aerosol loading at SMEAR Estonia lies between 0.08 and 0.21. As shown later, this relatively low ratio pertains to all the sites in this study: the maximum diffuse fraction of solar irradiance due to the direct aerosol effect was no more than 0.27 at the remote sites in boreal forests during the growing season.

Figure 4Ratio fdif_PAR∕fdif_bb as a function of fdif_bb (SMEAR Estonia). Different curves correspond to the best fits of the data for different solar zenith angles (sza).

Download

3.1.2 Parameterization of the diffuse fraction of PAR as a function of the diffuse fraction of broadband radiation

In this section we discuss the difference between the diffuse fractions of PAR and broadband radiation. Figure 4 displays the ratio between the diffuse fraction of PAR, fdif_PAR, and the diffuse fraction of broadband radiation,fdif_bb, as a function of fdif_bb where $fdif=R_{\mathrm{d}}/R_{\mathrm{g}}$. Since the radiation level depends on the solar zenith angle, we cast the daytime data into three solar zenith angle bins; the width of each bin was 10^∘. Each data set was fitted by the exponential function

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{fdi{f}_{\mathrm{PAR}}}{fdi{f}_{\mathrm{bb}}}}=a\exp (-(fdi{f}_{\mathrm{bb}}-b)/c)+d.\end{array}$$

The coefficients of the fitting function for each bin are reported in Table 2. Using this function, we can compare the diffuse fractions of PAR and broadband radiation over the whole range of sky conditions, including clear and cloudy skies. As expected (see more in Sect. 2.3), in Fig. 4 we observe an increase in the diffuse fraction of PAR up to 27 % compared with the value for broadband radiation at small fdif_bb, corresponding to clear-sky conditions. In absolute values, this difference between the diffuse fraction of the PAR and broadband radiation is not very large (e.g. fdif_bb=0.15 corresponds to fdif_PAR=0.18 under clear-sky conditions). However, as the diffuse fraction of global radiation in boreal forests varies in a relatively small range due to the direct effect of aerosol (e.g. Fig. 3), it is important to make corrections. As can be further noted from Fig. 4, the ratio fdif_PAR∕fdif_bb approaches one for overcast cloudy conditions, as in this case diffuse radiation prevails for both PAR and broadband radiation.

Table 2Best fit parameters for the diffuse fractions ratio, $\frac{fdi{f}_{\mathrm{PAR}}}{fdi{f}_{\mathrm{bb}}}$, as a function of the diffuse fraction of broadband radiation, fdif_bb, for various solar zenith angles (Eq. 7).

Download Print Version | Download XLSX

We use these results to obtain the diffuse fraction of global broadband radiation for the sites where only PAR was measured. In the following we use the term “diffuse fraction of global radiation” for broadband radiation.

Table 3Best fit parameters, correlation coefficients and p-values for radiation data ($R_{\mathrm{d}}/R_{\mathrm{g}}=k\mathrm{CS}+b$).

Download Print Version | Download XLSX

Figure 5The diffuse fraction of global irradiance as a function of CS (clear-sky data). The red symbols represent calculations with the clear-sky model, and blue symbols represent measurements.

Download

3.1.3 Aerosol influence on the diffuse fraction of global irradiance: comparative analysis for four sites

In this section, we consider the effect of aerosol on the diffuse fraction of global irradiance. In the following analysis the data were filtered to include only clear-sky conditions, based on the modelled and measured global irradiance, using Equation (7) for all of the sites. We deliberately used the equation based only on the global irradiance in order to demonstrate the effect of unfiltered cloud-contaminated data on the diffuse fraction of solar radiation. Figure 5 displays R_d∕R_g vs. CS at SMEAR I and II, SMEAR Estonia and Fonovaya (no aerosol data are available from Zotino). To separate the effect of clouds and aerosol particles, we report two quantities: the measured diffuse fraction of global irradiance (Fig. 5, blue symbols) and the modelled diffuse fraction of global irradiance (Fig. 5, orange symbols). Based on the analysis in Sect. 3.1.1, modelling provides information about the direct effect of aerosols on the diffuse fraction, while measurements illustrate the combined effects due to aerosols and clouds. Note that for consistency all the ratios R_d∕R_g were corrected in accordance with the previous section, meaning that only the ratios corresponding to the broadband radiation are reported (although PAR is measured at SMEAR I and II). For SMEAR I, the model data set includes 4 years, while only 2 years of measured data are available. Diffuse radiation is not measured at Fonovaya station; hence, only model results are shown. Moreover, the data from 2016 were not used due to forest fires in Siberia. Smoke plumes have large influences on the aerosol size distribution as a result of which the clear-sky model fails to predict diffuse and global radiation.

An increase in R_d∕R_g with increasing CS is observed at all of the sites, as follows from the model results (also representative for the measurements with an appropriate clear-sky equation, as discussed in Sect. 3.1.1 and demonstrated in Fig. 3). Note that the modelled values of R_d∕R_g correspond to the lower points in the measured data sets. The blue points above the modelled data are characterized by a larger diffuse radiation than those obtained for current AODs using Solis; hence, they represent the effect of clouds. According to the model calculations, the maximum diffuse fraction of global radiation due to the direct effect of aerosols did not exceed 0.27, while the minimum fraction was about 0.1. Furthermore, the aerosol population with CS smaller than approximately 0.005 s⁻¹ do not significantly contribute to light scattering (Fig. 5), as the diffuse fraction of global irradiance for these values of CS is almost constant and close to 0.1; this can be generally attributed to Rayleigh scattering.

We fitted modelled and measured data with the linear function $fdi{f}_{\mathrm{bb}}=k\mathrm{CS}+b$. The best-fit coefficients and correlation coefficients for four sites are reported in Table 3. All of the dependences pertaining to the modelled data, i.e. to the direct effect of aerosol particles on solar radiation, had correlation coefficients larger than 0.5 corresponding to a moderate correlation (except Fonovaya, where R=0.44, which was most likely due to the small data set). On the contrary, cloud-influenced data demonstrate rather weak correlations with $\mathrm{0.18}<R<\mathrm{0.33}$.

Figure 6Light use efficiency (LUE) as a function of the diffuse fraction of global irradiance (R_d∕R_g). All dependences are statistically significant (p<0.001).

Download

3.2 The effect of diffuse radiation on GPP

3.2.1 The effect of diffuse radiation on GPP: comparative analysis for all of the sites

In this section we study the LUE and PAR dependences on the diffuse fraction of global radiation in order to better understand the behaviour of the dependence of GPP on R_d∕R_g. Figure 6 displays the dependences of the LUE on R_d∕R_g for all sites. All of these dependences exhibit a linear relationship with the correlation coefficients between R=0.67 and R=0.83 (except Fonovaya with R=0.44, which can be attributed to both the short data set and the less precise gradient method used for the CO₂ flux calculations). The LUE slope reflects the canopy properties, i.e. it characterizes the ability of a forest stand to take up more CO₂ in response to an increasing diffuse fraction of solar irradiance. The steepest LUE slopes pertain to the mixed forests at SMEAR Estonia and Fonovaya, while the slopes are approximately 60 % less steep in coniferous forests (Table 4). This difference is presumably due to the forest type, as mixed forests have a larger potential for photosynthetic activity enhancement due to a larger leaf area index and a deeper canopy. We emphasize that the difference is seen in the LUE, in accordance with the LUE definition given in Sect. 2.5, which includes a dependence on LAI and tree height attributed to fAPAR in the standard definition. Note that the increase in the LUE from approximately 0.01 to 0.03–0.04 mol CO₂ mol photons⁻¹ observed for mixed forests is similar to that reported by Cheng et al. (2016) for mixed and broadleaf forests in the USA.

Table 4Linear regression coefficients for PAR and LUE at different sites: $\mathrm{LUE}=L_{\mathrm{1}}+L_{\mathrm{2}}\cdot (R_{\mathrm{d}}/R_{\mathrm{g}})$, $\mathrm{PAR}=R_{\mathrm{1}}+R_{\mathrm{2}}\cdot (R_{\mathrm{d}}/R_{\mathrm{g}})$.

Download Print Version | Download XLSX

Figure 7 displays the dependences of the global PAR on R_d∕R_g for all sites. As expected, PAR decreases as R_d∕R_g increases: at smaller values due to aerosol particles, and at larger values due to clouds. As follows from Fig. 5, the values of $R_{\mathrm{d}}/R_{\mathrm{g}}<(\mathrm{0.2}-\mathrm{0.27})$ mostly correspond to the influence of aerosol particles (but they can also be influenced by thin clouds), while larger values of R_d∕R_g are associated with the presence of clouds. Similarly to the LUE, these dependences are linear with high correlation coefficients ($\mathrm{0.78}<R<\mathrm{0.90}$). Generally, the slopes of the linear dependences in Fig. 7 were similar (within the range of 1081–1194 µmol s⁻¹ m⁻²), which can probably be attributed to similar cloud attenuating properties over all of the sites at middle latitudes. The exception is SMEAR I, where the slope is lower (944 µmol s⁻¹ m⁻²). Solar radiation under clear-sky conditions is also significantly lower at SMEAR I compared with the other sites, which is partly due to the high latitude, and partly because the growing season at SMEAR I is July and August (i.e. it does not include June which has the highest global irradiance values).

Figure 7Photosynthetically active global radiation (PAR) as a function of the diffuse fraction of global irradiance (R_d∕R_g). All dependences are statistically significant (p<0.001). The number of sample points is reported in the caption of Fig. 6.

Download

The vertical scattering of the data in Fig. 7 is presumably due to two factors: first, the variability in the radiation intensity during daytime and growing season, and second, the different influences of clouds, as the same diffuse fraction of global irradiance may pertain to the different attenuations of global radiation by clouds. Note that the latter factor is excluded from the Fonovaya data set by the parameterization. One can conclude that the PAR variability due to clouds was larger than the diurnal (associated with different solar zenith angles during the day) and day-to-day PAR variability in the growing season. Thus, additional binning by, e.g. solar zenith angle, would be redundant, as the decrease in PAR variability due to binning would be hidden by the stronger scattering due to clouds.

Finally, based on the linear dependences of the LUE on R_d∕R_g and PAR on R_d∕R_g, we can estimate how GPP depends on R_d∕R_g. When we multiplied the LUE by PAR, parabolic dependences were obtained for all the sites, with a maximum due to the effect of diffuse radiation on photosynthesis. Figure 8 shows the estimated GPP dependences on R_d∕R_g for the different sites for comparison, while Fig. 9 displays data sets and estimated curves separately for all sites, similar to Figs. 6 and 7.

Figure 8Estimated GPP dependences on R_d∕R_g for all of the sites (obtained as $\mathrm{GPP}=\mathrm{LUE}\cdot \mathrm{PAR}$ using the coefficients for PAR and the LUE dependences on R_d∕R_g reported in Table 4).

Download

Figure 9GPP dependences on R_d∕R_g for all of the sites. The curves represent estimated GPP (the same parabolas as in Fig. 8). We use a dashed curve for Fonovaya due to the relatively low correlation coefficient obtained for the LUE (R=0.44, Fig. 6). The data sets for all of sites were cast in bins in R_d∕R_g, the width of each bin is $R_{\mathrm{d}}/R_{\mathrm{g}}=\mathrm{0.04}$ ($R_{\mathrm{d}}/R_{\mathrm{g}}=\mathrm{0.08}$ for Fonovaya). Black points correspond to the mean GPP in each bin and error bars show the standard deviation of the data for each bin.

Download

3.3 Discussion

In this section we combine the results from the previous sections to make conclusions regarding the direct effect of aerosols on GPP, and we compare the results obtained with those from previous studies.

As previously mentioned in Sect. 3.1.3, a cloud-biased data set and a standard linear regression analysis results in weak but significant (p<0.001) correlations between CS and R_d∕R_g (Table 3). The relatively high cloud-biased diffuse fraction of global radiation at low CS leads to an underestimation of the effect of increasing aerosol loading for the cloud-biased data set. If CS increases from 0.002 to 0.015 s⁻¹ (obtained for the clear-sky conditions), a relative increase in the diffuse fraction of global radiation following from the clear-sky model is from 110 % to 165 % at all of the sites except Fonovaya, while this increase is between 65 % and 118 % for cloud-biased data. In the following we use only the results for the clear-sky model (representing the measured data set when the stricter Equation (8) of clear sky is applied, as follows from Fig. 3). In absolute values the increase was quite small: from 0.11 to ∼0.27 at SMEAR I and II, and from 0.11 to ∼0.2 at SMEAR Estonia and Fonovaya. The increases in R_d∕R_g over these value ranges led to increases in GPP from 17.2 to 18.6 µmol s⁻¹ m⁻² at Fonovaya, from 18.5 to 20.9 µmol s⁻¹ m⁻² at SMEAR Estonia, from 15.8 to 16.9 µmol s⁻¹ m⁻² at SMEAR II and from 8.8 to 10.0 µmol s⁻¹ m⁻² at SMEAR I. The largest relative increases in GPP due to the increasing aerosol loading from its minimum value to its maximum value were observed for SMEAR I and SMEAR Estonia (14 % and 13 % respectively). Note, however, that the median value of CS should be increased by a factor of about 5 at SMEAR I to get this maximum gain in GPP, whereas this same increase in GPP would be observed at SMEAR Estonia if the median CS increased by a factor of 2–3.

Overall, we obtained rather weak dependence of the diffuse fraction on CS. It is much weaker than that reported by Kulmala et al. (2014): for all of the sites the slope is less than 10 s (Table 3) compared with the almost 100 s obtained in the above-mentioned study. This difference is due to the inappropriate equation of clear sky selecting cloud-biased points with a diffuse fraction up to 0.8 (Kulmala et al., 2014) and a different statistical method (bivariate fitting compared with the linear regression used in this study). Note that Kulmala et al. (2014) reported minimum and maximum possible slopes for an increase in the diffuse fraction of global radiation with CS. Our present results are close to their minimum slope. Furthermore, due to the large diffuse fractions attributed to the effect of aerosols rather than clouds, the maximum direct effect of aerosols on GPP was overestimated by Kulmala et al. (2014). In the present study we obtained a 6 % increase in GPP at SMEAR II due to the diffuse radiation effect rather than the ≈30 % reported by Kulmala et al. (2014). However, their minimum slope, reported for GPP vs. R_d∕R_g dependence, would result in an increase of GPP similar to this study.

Note that the aerosol loading observed at all sites corresponds to 0.04 < AOD₆₇₅<0.35 with the typical values being in the range of 0.05–0.10 and AOD₅₀₀<0.25 (see Appendix A). In accordance with the study by Park et al. (2018), an increase in the diffuse fraction did not exceed 0.3 for these relatively low AOD values. Much higher diffuse fractions (0.5–0.7) due to the direct aerosol effect were obtained by Cirino et al. (2014) for the biomass burning season in the Amazon.

Next, all of the GPP dependences have a maximum due to clouds. The maximum corresponds to the clouds with a diffuse fraction in the order of 0.4–0.5. According to Cheng et al. (2016) and Pedruzo-Bagazgoitia et al. (2017), this R_d∕R_g corresponds to optically thin clouds with cloud optical thicknesses less than 5. Conversely, GPP decreases for optically thick clouds, which has also been demonstrated by Cheng et al. (2016). The largest increase is 32 %–33 % at SMEAR Estonia and Fonovaya, whereas the smallest increase is 11 % at SMEAR II compared with the GPP values on clear days characterized by low aerosol loading. At middle latitudes with a similar attenuation of radiation due to aerosols and clouds, the increase in GPP depends on the LUE slope: the steeper the LUE slope is, the more pronounced the maximum.

Based on Fig. 8, similar forest stands at Zotino and SMEAR I demonstrated a similar dependence of the GPP on R_d∕R_g, while this dependence was different for the coniferous forest at SMEAR II. GPP at SMEAR II under clear sky is almost 1.5 times larger than the corresponding GPP at SMEAR I and Zotino, but GPP increase under cloudy sky is smaller at SMEAR II. This could be a consequence of the closed canopy and the higher leaf area index of the SMEAR II forest stand. Our GPP data sets, reported in Fig. 9, look similar to those reported by Alton et al. (2007) and Alton (2008). The GPP dependence reported for SMEAR II is also similar to that reported by Alton (2008) for needle-leaf forests, but for mixed forests we obtained an increase of up to 30 % compared with the moderate 10 % increase for broadleaf forests reported by Alton (2008). Note that Alton (2008) used parameterization, Eq. (1), for the diffuse fraction of global radiation, while we had measurements of diffuse radiation at four sites out of five.

Finally, we considered the data from Zotino including the periods of forest fires (Park et al., 2018). Figure 7 suggests that forest fires do not have any specific influence on the PAR decrease with increasing R_d∕R_g compared with cloudy sky. In other words, plumes from forest fires lead to a similar decrease in PAR and a similar separation in diffuse and direct fractions as some clouds. The same holds for the LUE of an ecosystem: the dependence of the LUE on R_d∕R_g at Zotino is similar to that of other coniferous sites. However, a significant increase in GPP under wildfire plumes can potentially be obtained at a daily timescale because the radiation regime with (R_d∕R_g)_max, i.e. close to the optimal conditions for ecosystem photosynthesis, can persist for a long time under plume conditions. At the same time, clouds may be intermittent and the effect of a sporadic GPP increase can be compensated for by the smaller GPP when clouds are in front of the sun and radiation is reduced (Pedruzo-Bagazgoitia et al., 2017).