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

The effect of aerosol loading on solar radiation and the subsequent effect on photosynthesis is a relevant question for estimating climate feedback mechanisms. This effect is quantified in the present study using ground-based measurements from five remote sites in boreal and hemiboreal (coniferous and mixed) forests of Eurasia. The diffuse fraction of global radiation associated with the direct effect of aerosols, i.e. excluding the effect of clouds, increases with an increase in the aerosol loading. The increase in the diffuse fraction of global radiation from approximately 0.11 on days characterized by low aerosol loading to 0.2–0.27 on days with relatively high aerosol loading leads to an increase in gross primary production (GPP) between 6 % and 14 % at all sites. The largest increase in GPP (relative to days with low aerosol loading) is observed for two types of ecosystems: a coniferous forest at high latitudes and a mixed forest at the middle latitudes. For the former ecosystem the change in GPP due to the relatively large increase in the diffuse radiation is compensated for by the moderate increase in the light use efficiency. For the latter ecosystem, the increase in the diffuse radiation is smaller for the same aerosol loading, but the smaller change in GPP due to this relationship between radiation and aerosol loading is compensated for by the higher increase in the light use efficiency. The dependence of GPP on the diffuse fraction of solar radiation has a weakly pronounced maximum related to clouds.

1 Introduction

According to the IPCC report (Pachauri et al., 2014), the influence of aerosol loading on solar irradiance remains a relevant, open question. Aerosol particles influence the radiative balance of the Earth due to aerosol–radiation and aerosol–cloud interactions. The aerosol–radiation interactions include the scattering and absorption of solar radiation (direct effect), and potential changes in cloud properties associated with these radiative effects (semi-direct effect). Another important consequence of aerosol–radiation interactions is an increase in the diffuse fraction of solar radiation entering the Earth surface.

The influence of aerosol–radiation interactions on the diffuse fraction of solar radiation is relevant for estimates of the terrestrial carbon sink and for understanding climate feedback mechanisms (Kulmala et al., 2014). The increase in the land carbon sink due to enhanced industrial aerosols was estimated to be 20 %–25 % during the period of “global dimming” between 1960 and 1999 in the modelling study by Mercado et al. (2009). The physical mechanism behind the growth of the terrestrial carbon sink is as follows. An increased diffuse fraction of solar irradiance due to aerosols makes it easier for light photons to penetrate into the canopy. Moreover, diffuse light coming from different angles increases the efficiency of CO₂ uptake by leaves of different orientation (Alton et al., 2007). This leads to an increase in the light use efficiency (LUE) of plants, which quantifies the amount of CO₂ fixed by an ecosystem per unit of photosynthetically active radiation (PAR), and gross primary production (GPP), which quantifies the amount of CO₂ fixed by an ecosystem per unit area per unit time. However, this mechanism is ecosystem-dependent and presumably depends on canopy height and the leaf area index (LAI) (Cheng et al., 2015; Kanniah et al., 2012; Niyogi et al., 2004): an enhanced diffuse radiation does not lead to an increase in grasslands' GPP, while GPP of broadleaf forests can be significantly increased. Similarly, a study based on AmeriFlux data from several sites (Cheng et al., 2015), including broadleaf forests and crops, suggested an increase in GPP due to diffuse radiation for forests but not for crops.

Estimates of the aerosol effect on GPP are uncertain for two reasons. First, it is not clear how large the effect of aerosols is on diffuse radiation. Second, the associated effect of diffuse radiation on GPP can be both negative and positive (Alton, 2008; Niyogi et al., 2004; Park et al., 2018). Therefore, aerosol–radiation interactions may either lead to an increase or a decrease in GPP, depending on aerosol loading. As an example, Niyogi et al. (2004) reports an increase in GPP for broadleaf forests at any aerosol loading typical for the sites they considered. On the contrary, a recent model study suggests that for a substantial proportion of boreal forests in the territories of Finland, Estonia and Russia, the direct aerosol effect at relatively high aerosol loading leads to a decrease in the annual diffuse irradiance and GPP (Lu et al., 2017). Estimates on the effects of solar radiation on GPP have often been obtained by utilizing parameterizations (Alton, 2008) or based on the results of numerical modelling and satellite observations (Lu et al., 2017). The aims of this study are to provide a comprehensive data analysis related to the direct effect of aerosol on solar radiation, to separate cloud and aerosol effects on solar irradiance, and to further quantify the influence of solar radiation on GPP. The data sets include continuous field measurements from five stations in Finland, Estonia and Russia.

Note that this analysis can be put into a broader context regarding the quantification of terrestrial feedback loops, e.g. the COntinental Biosphere–Atmosphere–Cloud–Climate (COBACC) feedback loop (Kulmala et al., 2014). This loop considers the effect of the carbon dioxide concentration on biogenic volatile organic compound emissions (BVOC), further relates BVOC-aerosol interactions to the variability in solar radiation and finally the loop is closed with the effects of radiation for the ecosystem GPP and CO₂ concentration. Here we provide a quantification of the part of this loop related to aerosol–solar radiation–photosynthesis interactions in boreal and hemiboreal forests.

In line with the aims of this study, we focus on two problems: (1) aerosol–radiation interactions, for which we quantify the diffuse fraction of solar radiation that can be observed due to the direct aerosol effect, and (2) the diffuse radiation effect on photosynthesis. First, we quantify the effect of aerosol–radiation interactions on diffuse radiation in boreal forests. For the estimates of the effect of aerosol on solar radiation it is important to separate clear times of the day from cloudy periods, because clouds are much more effective than aerosols with respect to scattering solar irradiance. In some previous studies (Gu et al., 2002; Kulmala et al., 2014) this separation was made based on the ratio of measured global irradiance to theoretically calculated global irradiance at the top of the atmosphere. This criterion of clear days can indeed be acceptable when one needs to distinguish between mostly sunny/mostly cloudy days with respect to the incoming global irradiance or the incoming solar energy. However, it fails to separate sunny/cloudy periods from the point of view of the amount of diffuse irradiance, which is crucially important for photosynthesis. Generally, radiation data used for analyses is averaged over a time interval (often 30 min). Some types of clouds, such as cumulus, are rather intermittent with oscillation periods of several minutes (Duchon and O'Malley, 1999). During the periods when there is no cloud in front of the sun, global irradiance can be even higher than theoretically predicted due to a global radiation enhancement (Pecenak et al., 2016). As a result, an averaged global irradiance can be close to that theoretically predicted for the clear sky, while averaged diffuse irradiance is significantly enhanced at the same time due to the presence of clouds. Such data can be erroneously attributed to clear-sky conditions, meaning that the effects of clouds on the diffuse radiation can be associated with the direct aerosol effect. These kinds of conditions are most likely to cause the highest GPP, as global irradiance does not decrease while diffuse irradiance is already high. Second, we consider the effect of solar radiation on GPP. We investigate the effect of the diffuse fraction of solar radiation on the LUE and further quantify the maximum effect on GPP due to aerosol–radiation interactions for different ecosystems in a boreal forest.

2 Materials and methods

In this section we introduce data sets and methods used in the study. In Sect. 2.1 we present the sites and data sets. Section 2.2 describes the clear-sky model used to separate clear sky and clouds. We discuss the consequences of having different (PAR or broadband) radiation measurements at different sites in Sect. 2.3 and suggest a method to make the data sets comparable. In Sect. 2.4 we introduce the condensation sink as a measure of aerosol loading. Finally, Sect. 2.5 describes the method used to study the diffuse radiation effect on ecosystem GPP using a separation of GPP into the LUE and the PAR.

Table 1 Data sets from different sites used in this study.

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 1 Location of the sites (see Table 1 for latitudes and longitudes of the different stations).

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

$$\begin{array}{}\text{(1)}& fdi{f}_{\mathrm{bb}}=R_{\mathrm{d}}/R_{\mathrm{g}}=\mathrm{1.45}-\mathrm{1.81}x,\end{array}$$

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:

$$\begin{array}{}\text{(2)}& \mathrm{GPP}=\mathrm{TER}-\mathrm{NEE},\end{array}$$

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:

$$\begin{array}{}\text{(3)}& I_{\mathrm{gh}}=R_{g,\mathrm{mod}}=I_{\mathrm{0}}^{\prime }\exp \left(-{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{g}}}{{\cos }^{\mathrm{g}}\left(sza\right)}}\right)\cos \left(sza\right),\end{array}$$

while the direct irradiance, I_dir, was calculated as

$$\begin{array}{}\text{(4)}& I_{\mathrm{dir}}=I_{\mathrm{0}}^{\prime }\exp \left(-{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}}}{{\cos }^{\mathrm{b}}\left(sza\right)}}\right).\end{array}$$

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:

$$\begin{array}{}\text{(5)}& I_{\mathrm{dh}}=R_{d,\mathrm{mod}}=I_{\mathrm{gh}}-I_{\mathrm{dir}}\cos \left(sza\right).\end{array}$$

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):

$$\begin{array}{}\text{(6)}& \mathrm{CS}=\mathrm{2}\mathit{\pi }{D}_{\mathrm{v}}\underset{\mathrm{0}}{\overset{d_{\mathrm{p},\max }}{\int }}{d}_{\mathrm{p}}\mathit{\beta }n\left(d_{\mathrm{p}}\right)d{d}_{\mathrm{p}},\end{array}$$

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 Results and discussion

We present the results of our study in two subsections. In Sect. 3.1 we report the results related to the aerosol effect on solar radiation, and in Sect. 3.2 we report the results related to the effect of diffuse radiation for ecosystem photosynthesis. The link between these results and their relation to other studies are discussed in Sect. 3.3.

Figure 2 Modelled vs. measured irradiance for SMEAR Estonia: measured global radiation, R_g,meas; modelled global radiation, R_g,mod; measured diffuse radiation, R_d,meas; and modelled diffuse radiation, R_d,mod. (a) 1 June 2016 – clear day; (b) 6 June 2016 – cloudy day. Timescale corresponds to local winter time (UTC+2).

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 3 The 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.

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 4 Ratio 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).

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 2 Best 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).

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 3 Best fit parameters, correlation coefficients and p-values for radiation data ($R_{\mathrm{d}}/R_{\mathrm{g}}=k\mathrm{CS}+b$).

Figure 5 The 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.

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 6 Light 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).

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 4 Linear 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}})$.

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 7 Photosynthetically 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.

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 8 Estimated 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).

Figure 9 GPP 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.

3.2.2 Constraints on the LUE and the diffuse fraction of solar radiation associated with the maximum ecosystem GPP under diffuse light

Well-pronounced linear dependences of the LUE and PAR on R_d∕R_g can be used to estimate how large an increase in the LUE should be in order to have GPP increase under diffuse radiation and at what diffuse fraction of solar radiation the maximum GPP can be observed. If

$$\begin{array}{}\text{(10)}& \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}}),\end{array}$$

then the maximum GPP is reached at $(R_{\mathrm{d}}/R_{\mathrm{g}}{)}_{\max }=-\mathrm{0.5}(L_{\mathrm{1}}/L_{\mathrm{2}}+R_{\mathrm{1}}/R_{\mathrm{2}})$, estimated as the point where the parabola $\mathrm{GPP}=\mathrm{LUE}\cdot \mathrm{PAR}$ has its maximum. The position of this maximum depends on the ratios L₂∕L₁ and R₂∕R₁. For a certain range of parameters, the maximum of the parabola can be located at $R_{\mathrm{d}}/R_{\mathrm{g}}<\mathrm{0.08}$, which is below the minimum diffuse fraction measured at our sites; therefore, it is not feasible for the latitudes we consider in this study. In this case, GPP monotonically decreases when R_d∕R_g increases from ∼0.1 to 1. Conversely, (R_d∕R_g)_max should be larger than 0.08–0.1 for the GPP to have a maximum under diffuse light. Note, that PAR dependences on R_d∕R_g at middle latitudes are similar: $R_{\mathrm{1}}/R_{\mathrm{2}}\approx -\mathrm{1.5}$, while for SMEAR I this ratio is $R_{\mathrm{1}}/R_{\mathrm{2}}\approx -\mathrm{1.3}$. From these estimates, $L_{\mathrm{2}}>L_{\mathrm{1}}/\mathrm{1.2}$ for the middle latitudes. Since L₁ is roughly the minimum value of the LUE at $R_{\mathrm{d}}/R_{\mathrm{g}}\approx \mathrm{0.1}$ (clear sky), while L₁+L₂ is the maximum of the LUE at $R_{\mathrm{d}}/R_{\mathrm{g}}=\mathrm{1}$ (overcast conditions), L₂ can be treated as the maximal gain in the LUE under diffuse light. Thus, the GPP will have a maximum associated with diffuse radiation if the ecosystem LUE under diffuse light increases by more than approximately 80 % of its minimum possible value (which is observed under clean conditions on clear days). For the sites considered in the present study, the smallest gain in the LUE due to diffuse radiation is observed at SMEAR II, where the LUE under diffuse light was almost twice as large as its value on clear days. The largest gain was at SMEAR Estonia and Fonovaya where the LUE grew by almost a factor of 3 if the dominating radiation conditions in the area changed from mostly direct to mostly diffuse radiation. Therefore, all ecosystems displayed maxima of GPP dependence on R_d∕R_g due to diffuse light, although at different values of R_d∕R_g.

Moreover, this approach clearly demonstrates that the maximum GPP can never be reached under overcast conditions. If we again take $R_{\mathrm{1}}/R_{\mathrm{2}}=-\mathrm{1.4}$ as for the middle latitudes, then the position of the maximum is at $(R_{\mathrm{d}}/R_{\mathrm{g}}{)}_{\max }=-\mathrm{1}/\mathrm{2}(L_{\mathrm{1}}/L_{\mathrm{2}})+\mathrm{0.7}$. One can immediately deduce that for large slopes of the LUE, i.e. when L₁∕L₂ approaches zero, (R_d∕R_g)_max approaches 0.7. At SMEAR I, this position is restricted by (R_d∕R_g)_max ≈ 0.65. The maximum of GPP parabolas for the five sites considered in this study is at $(R_{\mathrm{d}}/R_{\mathrm{g}}{)}_{\max }\approx \mathrm{0.4}$–0.5.

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).

4 Conclusions

We quantified the direct effect of aerosol on solar radiation and GPP in boreal and hemiboreal forests in Eurasia. The analysis was based on data from five sites including coniferous and mixed forest ecosystems. The diffuse fraction of global radiation due to the direct aerosol effect was estimated to be in the range of $\mathrm{0.11}<R_{\mathrm{d}}/R_{\mathrm{g}}<\mathrm{0.27}$ at all of the sites.

For the first time we demonstrated a connection between solar radiation properties (the diffuse fraction of global radiation) and condensation sink. The latter parameter is used in aerosol studies and it is obtained from ground-based observations. Employing CS instead of a column-averaged aerosol parameter AOD is a necessary step towards further investigation of the COBACC climate feedback loop, linking biogenic volatile organic compounds emissions and aerosol characteristics.

The GPP-radiation analysis was performed using the separation of GPP into the LUE and PAR. We found a linear dependence between the diffuse fraction of solar radiation and the LUE, as well as between the diffuse fraction of solar radiation and PAR, for all of the sites. While the PAR dependences were quite similar to one another (except for SMEAR I which is located at a relatively high latitude), the LUE dependences were different: the slopes were 60 % steeper for mixed forests than for coniferous forests, and the intercepts were about 40 % lower for coniferous forests with open canopies. We obtained a parabolic shape for the GPP dependence on the diffuse fraction of solar radiation. The maximum of the parabola was more pronounced for mixed forests due to the above-mentioned differences in the LUE dependences between the mixed and coniferous forests. Note that parabolic, or near parabolic, shapes have been reported for different forest sites by Alton (2008) and Moffat et al. (2010) using different methods to those used in this study.

We showed that GPP can be increased by 6 %–14 % due to the direct effect of aerosol particles at remote sites compared to clean conditions with low values of CS. The maximum increase was observed for mixed forests at mid-latitudes and for coniferous forests at relatively high latitudes.

Furthermore, based on the similarity in the PAR dependences on the diffuse fraction of solar radiation for all of the sites, we obtained the constraints on the ecosystems' LUE increase under diffuse light necessary for a GPP maximum due to diffuse light. At mid-latitudes, the LUE of an ecosystem should increase by more than ∼80 % under diffuse light compared with its value under clear-sky conditions. Moreover, at the mid-latitude sites, the diffuse fraction of solar radiation corresponding to the maximum GPP can not exceed 0.7.

The specific shape of the GPP dependence on the diffuse fraction of solar radiation suggests that clouds with a 30 min-averaged fraction R_d∕R_g between 0.4 and 0.5 play an important role in determining ecosystems' GPP and demand further investigation. An increase in GPP due to clouds can reach 32 %–33 % for mixed forests and 21 %–26 % for coniferous forests with an open canopy. Other relevant questions include cloud effects on the radiation regime and ecosystems' GPP on annual scale and the investigation of potential aerosol effects on the evolution of clouds over forests.

Appendix A: Condensation sink vs. AOD

To make our study comparable to studies that use AOD₅₀₀ as a parameter quantifying aerosol loading, we performed an additional analysis for the data sets from SMEAR II and SMEAR Estonia. We used 30 min-averaged AOD from AERONET sites and CS, corresponding to the period from 09:00 to 15:00 (local time) and maximal growing season. The data set from SMEAR II encompasses 3 years (2008–2010) and the data set from SMEAR Estonia encompasses 4 years (2013–2016).

The results are shown in Fig. A1. Though AOD₅₀₀ clearly increases with increasing CS, the scatter of data is great. It means that in spite of the fact that both parameters are roughly proportional to the aerosol surface area distribution (Sundström et al., 2015), and in spite of the presumably well-mixed boundary layer during daytime, there is no simple relationship between ground-based and column-integrated aerosol characteristics. This has also been noted by Sundström et al. (2015) for sites in South Africa. Sundström et al. (2015) found a moderate correlation between the in situ measured scattering coefficient and AOD from AERONET. At the same time, the scattering coefficient was strongly correlated with CS; therefore, a moderate correlation between CS and AOD from AERONET could be expected. Figure A1 shows a moderate correlation (R=0.53) for boreal forests at the mid-latitudes.

Figure A1 Aerosol optical depth at 500 nm (AOD₅₀₀) as a function of the condensation sink (CS) for SMEAR II and SMEAR Estonia (R=0.53, p<0.001).