Impact of Atmospheric Optical Properties on Net Ecosystem Productivity of Peatland in Poland

Peatlands play an important role in the global carbon cycle due to the high carbon storage in the substrate. Ecosystem production depends, for example, on the solar energy amount that reaches the vegetation, however the diffuse component of this flux can substantially increase ecosystem net productivity. This phenomenon is observed in different ecosystems, but the study of the atmosphere optical properties on peatland production is lacking. In this paper, the presented methodology allowed us to disentangle the diffuse radiation impact on the net ecosystem production (NEP) of Rzecin peatland, Poland. It allowed us to assess the impact of the atmospheric scattering process determined by the aerosol presence in the air mass. An application of atmospheric radiation transfer (ART) and ecosystem production (EP) models showed that the increase of aerosol optical thickness from 0.09 to 0.17 caused NEP to rise by 3.4–5.7%. An increase of the diffusion index (DI) by 0.1 resulted in an NEP increase of 6.1–42.3%, while a DI decrease of 0.1 determined an NEP reduction of −49.0 to −10.5%. These results show that low peatland vegetation responds to changes in light scattering. This phenomenon should be taken into account when calculating the global CO2 uptake estimation of such ecosystems.


Introduction
The interaction between terrestrial ecosystems and the Earth's atmosphere is one of the most critical ecological issues, since the properties of the atmosphere determine each ecosystem's sustainability, development, or extinction [1]. The water and carbon dioxide balances of the biosphere are determined by both biotic and abiotic factors, and their global cycles are efficiently controlled by the solar energy that flows throughout the Earth's system. Thus, physical parameters of the atmosphere such as temperature and transmissivity determine the global carbon cycle. Changes in the physical properties of the atmosphere determine the fate of the ecosystems, and their carbon balance are important issues in the context of the greenhouse effect phenomenon. Earth surface-atmosphere CO 2 exchange is strongly related to photosynthesis and respiration rates. Both processes depend directly and indirectly on the total (SWin) solar radiation (the sum of direct (SWdir) and diffuse (SWdif) radiation) [2]. Thus, research into the optical parameters of the atmosphere is a basis for the assessment of one of the most fundamental controls of CO 2 exchange.
In the short time scale, the ecosystem's C balance consists of the net ecosystem production (NEP) which is the result of two key processes: gross primary production (GPP) and ecosystem respiration (ER) (Equation (1)) [3]: where GPP represents the vegetation's CO 2 uptake through photosynthesis, while ER is the release of CO 2 into the atmosphere through autotrophic and heterotrophic respiration processes.
Since peatlands play one of the most important roles in the global carbon cycle, as an ecosystem with significant potential for a highly positive impact on the climate system (net cooling effect over centuries) [4][5][6], the study of the CO 2 exchange between these ecosystems and the atmosphere is crucial for the prediction of the functioning of this huge natural soil carbon stock under global warming that appears in all climate zones [7][8][9][10]. This C storage has been formed through the photosynthetic fixation of carbon dioxide by vegetation, which in the long term has been larger than the release of C through plant respiration and peat decomposition [11]. Thus, the transformation of peatlands into carbon sources as the result of global warming is currently one of the biggest threats that may contribute substantially to the instability of the Earth-atmosphere thermodynamic system. Simultaneously, these ecosystems may prevent the loss of stability caused by crossing the planetary threshold. In other words, they have the potential to stop the independent process of increasing the greenhouse gas concentration and significantly slowing down or even reversing these couplings [12].
The canopy net photosynthesis is found to be dependent not only on the quantity but also on the degree of scattering of solar radiation, which has been confirmed by numerous simulations [13][14][15]. The increased production under more diffuse radiation is explained inter alia, as diffuse radiation more effectively penetrates the plant canopy [16][17][18]. This improves the photosynthetic activity of the lower located leaves in the canopy [2]. Longterm observations show that in the last several decades, both positive (dimming -until mid-1990s) and negative (brightening) trends of diffuse irradiance have been observed on the Earth's surface [19][20][21][22]. The trends of total solar radiation have been found as well [21], and their significance and direction are often related to local circumstances, e.g., urbanization [23][24][25]. Light scattering processes in the atmosphere are mainly determined by the presence of suspended components in the air, e.g., gases and different sized particles in the solid and liquid phases [26]. The presence of atmospheric aerosols, both at their height of occurrence in the air column [27] and the aerosol type [28,29], can affect the climate system due to their impact on the heat and water balances at both local and global scales [30]. Since aerosols also play an important role in the scattering process of the atmosphere [31,32], atmospheric radiative transfer (ART) models have been developed (e.g., [33]) in order to effectively simulate the radiation processes in the atmosphere [34]. The calculation of the solar irradiance with this model is realized on the basis of the parameters that integrate the effect of the presence of the aerosol on air properties. The aerosol optical thickness (AOT) is defined as the parameter that represents the airborne aerosol loading in the atmospheric column [35]. This measure can be effectively used for the assessment of the radiation transfer processes in the atmosphere [36]. AOT values measured on the ground are usually realized with sun photometers, which are often deployed within networks that provide technical support in terms of calibration and data processing, e.g., aerosol robotic networks (AERONET) [37].
The impact of diffuse solar radiation on ecosystem production was initially analyzed for agricultural areas, e.g., soybean, maize, peanut, and sugar beet crops [38][39][40][41]. Later studies were extended to other types of ecosystems, mainly deciduous and coniferous forests, tundra, grassland, and the Canadian peatlands [14,16,[42][43][44][45][46][47][48][49][50], but there is still a lack of such analyses for central European peatlands. The scattering process that determines the diffuse radiation amount is the result of the presence of both clouds and aerosols in the atmosphere. However, in a recently published study of diffuse radiation impact on peatland productivity, the separate effects of cloud and aerosol presence in the air column were not presented. In the case of the Canadian peatland studies, the analyses were focused on the impact of cloudiness on ecosystem production using a clearness index [50] or light use efficiency models [15,51], while New Zealand bog studies were based on different light response functions [52]. Thus, the main goal of the study was a quantitative estimation of the impact of optical parameter changes, determined by the aerosol presence in the atmosphere, on peatland net CO 2 uptake in cloudless conditions.

Site Description
Data presented in this paper were obtained at Rzecin peatland, located in the northwestern part of the Greater Poland Region, Poland, in Rzecin village (52 • 45 N; 16 • 18 E; 54 m a.s.l) (Figure 1) [53]. The peatland is located around 70 km northwest of Poznań city in the Notecka Primeval Forest complex. This peatland is classified as fen and the vegetation is comprised of typical peatland species, e.g., roundleaf sundew (Drosera rotundifolia), bog cranberry (Oxycoccus palustris), and bog sedge (Carex limosa). There are also various mosses found at this site with rigid bog moss (Sphagnum teres) (Schimp.) being the dominant species [54]. Tall vascular plants, e.g., common reed (Phragmites australis), broadleaf cattail (Typha latifolia), and sedges (Carex spp.) with LAI >4.8 m 2 ·m −2 dominate at the edge of the peatland, while bog mosses (Sphagnum spp.), fine bogmoss (Sphagnum angustifolium), flat topped bog moss (S. fallax), rigid bog moss (S. teres), and short, vascular plants like bog cranberry (Oxycoccus palustris), bog sedge (Carex limosa), bottle sedge (C. rostrate), wollyfruit sedge (C. lasiocarpa), common cottongrass (Eriophorum angustifolium), and roundleaf sundew (Drosera rotundifolia), with LAI < 1.0 m 2 ·m −2 dominate in the middle of the ecosystem [55,56]. At the eastern part of this ecosystem, there is a shallow lake (remnants of a larger body of water) that has been partly covered within the peat formation process [57]. The floating carpet of poorly degraded peat (approximately 70 cm thick) is a particular feature of the Rzecin peatland. This is an effect of the moss layer growing process that has been ongoing on the lake surface over the last 200 years [58]. There are cultivated fields, pastures, and pine forests located around the peatland, and their impact on the ecosystem has been noticeable over the last 200 years [59]. The site substrate water pH is 4.9, and conductivity is 52.8 S·cm −1 [60]. A Nature 2000 site was established at the site in 2006 due to both plant and bird species richness. The site is characterized by a moderately warm climate, where the annual mean air temperature is 8.5 • C, the annual mean precipitation is 526 mm, and there is a prevailing western wind [61].

Meteorological Conditions
Between 13 June and 18 September, the average 30-minute air temperature was 17.6 • C (minimum 2.0 • C, maximum 35.3 • C), while the vapor pressure deficit ranged from 0 to 36 hPa, and the average relative humidity over this time was 79%-with minimum and maximum values of 32% and 100%, respectively. The DI values varied between 0.12 and 1.0 over the studied period. The mean value of AOT500 collected with CIMEL sun photometer was equal to 0.15234, while the minimum and maximum values were 0.03535 and 0.56571, respectively. The sum of total solar radiation and precipitation for this period were 1745 MJ·m −2 and 286.6 mm, respectively.

Eddy Covariance
The eddy covariance technique (EC) is a widely used method of CO 2 exchange estimation [2,62] on the scale of the entire ecosystem (100-2000 meters) [63,64] without vegetation and soil disturbances [62,65]. This method is based on the simultaneous measurement of instantaneous values (fluctuations) of mass, e.g., CO 2 concentration and vertical wind speed component. Such a measurement strategy allows the direct observation of vertical mass exchange between the atmosphere and ecosystem surface [66]. EC technique is currently the state-of-the-art measuring method that provides the most reliable estimation of gas fluxes under field conditions [62].
The measurements of CO 2 exchange at the Rzecin site were carried out with the eddy covariance system installed on a 4.5 m tall tower that was located in the middle of the peatland [53]. The EC system consists of a sonic anemometer R3-100 (Gill Ltd., UK) and an enclosed CO 2 /H 2 O gas analyzer LI-7200 (LI-COR, USA), which measures the fluctuations of the vertical wind speed component and CO 2 /H 2 O concentration, respectively, at 20 Hz frequency.
The 30-minute values of the CO 2 fluxes were calculated with EddyPro [67]. Initially, the quality of the raw time series statistical tests recommended by [68] were applied (amplitude resolution, dropout, absolute limit, skewness, and kurtosis), as well as the spike count/removal according to [69]. Then, the following procedures were applied: compensation of the fluctuation of air density [70], correction for spectral losses [71], flagging according to Carbo Europe standard [72], and footprint calculations [73]. The footprint analysis showed that the origin of the measured fluxes was located within a 200 m radius around the tower.

Sun Radiometer
The aerosol optical thickness (AOT) was measured by the multiband automatic sun/sky radiometer CE318 (CIMEL Electronique, Paris, France) that has been operating at the study site since May 2016. This solar-powered and weather-resistant robotic device has 9 spectral bands [37,74] and it is installed at a dedicated tripod tower (52 • 45 43.2 N; 16 • 18 34.4 E) ( Figure 1) at a height of 4 m above the surface a short distance from EC tower. The data obtained with the CE318 is automatically processed and transmitted to the AERONET website [75]. This ground-based global monitoring system coordinated by NASA consists of around 900 measurement stations around the world [74][75][76]. The strength of this network relies on provision and standardization of the measurement protocol, data processing, and calibration, allowing for multiyear and large-scale comparisons [37].

Meteorological Measurements
The measurements of the air temperature (Ta) and relative humidity were carried out using a HC2A-S humidity probe installed at a standard height of 2 m above the surface (RotronicMessgeräte GmbH, Ettlingen, Germany). The values of both total (I t ) and diffuse (I f ) photosynthetic photon flux density (PPFD) were obtained with the sunshine sensor BF5 (Delta-T Devices, Burwell, UK) [77]. It was installed at a height of 3.0 m, and was placed on the same scaffold as the EC system. The height of BF5 installation prevented the obstacles e.g., tall plants in the sunshine sensor field of view. The diffusion index (DI) was calculated on the basis of these values as the ratio of I f and I t [47]. The solar radiation was measured also with a pair of upward and downward facing CM3 pyranometers (Kipp&Zonen, Delft, the Netherlands) that were used for total (SWin) and reflected (SWref) shortwave radiation flux densities, respectively. An additional pair of Quantum sensors (SKP215, Skye Instruments Ltd., Powys, UK) was applied for the measurements of both total (I t ) and reflected (I ref ) photosynthetic photon flux density. Both pairs of sensors were installed on a steel arm at a height of 2.35 m above the Menyantho trifoliatae-Sphagnetum teretis Warén 1926 vegetation community [78], next to EC tower. This set of sensors provided the data that was used for the Broad-band Normalized Difference Vegetation Index (NDVIb) calculation according to the following formula (Equations (2)-(4)) [79]: where: and: , RSin is total shortwave radiation flux density [W·m −2 ], and RSref is reflected shortwave radiation flux density [W·m −2 ]. The values of NDVIb correlate well with satellite normalized difference vegetation index (NDVI) data, but it does not suffer from discontinuity problems that stem from the remotely obtained observations [80]. Thus, the NDVIb was used here as a proxy for vegetation CO 2 uptake capacity since the relation between NDVIb and gross ecosystem production (GEP) was found to be noticeable on the studied site [81]. Within this study, the gross ecosystem production (GEP) is approximated as GPP, and small discrepancies between GPP and GEP were ignored here.

Ecosystem Productivity (EP) and Atmospheric Radiative Transfer (ART) Modeling
There are several ecosystem production models developed for peatlands [82][83][84][85]. They are effective tools for CO 2 uptake estimation, but the solar radiation scattering process is not taken into account as a factor in these algorithms. There were several models developed for higher vegetation (usually trees) where the spatial structure is taken as a factor that impacts the solar radiation interception in the plant canopy [86,87]. Due to the limited information about the peatland canopy structure, the simplified diffuse radiation ecosystem production (EP) model was applied within this study, where NEP was estimated on the basis of the separate impact of direct (I r ) and diffuse (I f ) solar radiation [87] (Equation (5)). This model consisted of two ecosystem CO 2 balance components: GEP, where the rectangular hyperbola equation's parameters α (initial quantum yield) and β (empirical coefficient) are estimated separately for I r and I f , while ER is estimated on the basis of air temperature (Ta): where NEP is net ecosystem production flux [µmol·m −2 ·s −1 ], Ta is air temperature at 2 m height [ • C], I f is diffuse PPFD [µmol·m −2 ·s −1 ], I r is direct PPFD [µmol·m −2 ·s −1 ], I t is total PPFD [µmol·m −2 ·s −1 ], and c 1 , c 2 , α f , α r , β f , and β r are empirical parameters.
The ecosystem respiration flux density (the exponential part of the formula) is assumed to be Ta dependent. A simple atmospheric radiative transfer (ART) model described by [33] has been applied to simulate solar irradiation ( Figure 2). This model is mostly used for renewable energy applications [88], for ecosystem productivity simulations [89], and for cloud screening data of solar irradiation [90]. Some empirical formulas described by Justus and Paris, [33] are used to define other ATM models [91,92]. Although the Justus and Paris ART model is very simple and it shows relatively good agreement with clear-sky observations (root mean square error of 3.5%). In this study, we used the second class pyranometer (CM3) with an accuracy of ±10% for daily sum flux, therefore the use of the simple ART model is justified. This ART model calculates clear-sky spectral direct and diffuse irradiance, spectral absorption within the atmosphere, and the upward reflected spectral irradiance at the top of the atmosphere. The irradiance model, based on similar approaches by [93], evaluates the spectral irradiance between 0.29 and 4.0 µm, with a resolution that varies from 0.005 to 0.1 µm. This empirical ART includes absorption by two uniformly mixed trace gases (water vapor and ozone) as well as both the scattering and absorption by aerosol and the scattering by molecular particles. Both molecular and aerosol optical properties are modeled with a simple wavelength-dependent optical thickness described by the Angstrom exponent, single scattering albedo, and asymmetry parameter. In the case of ozone and water vapor, the total columnar values are assumed. The broadband surface albedo is used to estimate the multi-scattering effect of incoming and outgoing solar radiation. There are the following input parameters used for testing and parameterizing this model: Julian day, solar zenith angle, aerosol optical thickness at 500 nm, total water vapor [g·cm −3 ], single scatter albedo at 500 nm (SSA500 = 0.95), asymmetry parameter (g = 0.65), Angstrom exponent, and surface albedo. The values of AOT500 were applied during all calculations because this wavelength is in the range of radiation used for the photosynthesis process.
The compilation of ART and EP models was performed in order to develop the model where NEP is determined by the AOT values. In this case, the output values of direct and diffuse irradiance from the ART model are used as input data for the EP model ( Figure 2).

Broad-band Normalized Difference Vegetation Index
The first step of data selection consisted of the extraction of the period when the plant's photosynthetic efficiency was stationary in order to avoid the canopy stage of development impact on CO 2 exchange parametrization. Thus, the NDVIb's weekly populations collected during 2016 were analyzed with a post hoc Tukey's test to identify the periods of similar NDVIb values.

Vapor Pressure Deficit (VPD)
Due to the fact that the daytime leaf stomata activity is vapor pressure deficit dependent [94], its impact on NEP was analyzed within the second step of data selection. It was found that if VPD value exceeds a certain threshold, it determines the stomata closing and the gas exchange between the plant and ambient air is substantially reduced or stopped [95,96].

Clear Sky Conditions
In order to estimate the separate aerosol effect on the peatland productivity, the clear sky conditions were needed to be extracted from the previously selected data set. Since the AOT values collected by the sun photometer are processed using the quality assessment and quality check (QAQC) algorithm by AERONET, the final AOT database contains three levels of quality data: 1.0, 1.5, and 2.0. All unprocessed data obtained by the radiometer are stored on level 1.0, while level 1.5 is selected by the cloud-screening and quality control algorithm application [37,75]. The data from level 2.0 is the result of corrections that are based on AERONET calibration procedures and reflect cloudless conditions. The periods when AOT data from level 2.0 were available were used as the last criteria of the data selection. Finally, the selected set of NEP, AOT, and meteorological data was the basis of parametrizing EP and ART models in order to estimate the Rzecin peatland CO 2 productivity change to the AOT modifications during clear sky conditions. The simulation was performed on 23 June 2016 since it was characterized by clear sky conditions, all necessary data availability, optimal plant development stage, and VPD values below 20 hPa. Additionally, it is one of the longest days of the year.

EP and ART Models Accuracy Assessment
The models were evaluated in terms of accuracy using the following parameters: mean absolute error (MAE), root mean square error (RMSE) [97], and normalized root mean square error (NRMSE). The standard error (SE) of EP model parameters was calculated according to the following formula (Equation (6)) [98]: where σ is population standard deviation and N is the number of elements in the sample. All analyses presented in this paper were performed within R software [99].

Seasonal Patterns of NDVIb
The NDVIb populations were found to be homogeneous within the period from week 24 to week 39 of 2016 ( Figure 3). The test indicated no statistical difference (significance level p = 0.05) between the NDVIb population of week 24 and the populations placed within the weeks from 25 to 39. On the basis of this analysis, the assumption that the CO 2 uptake capacity of the plant canopy was constant over this period was done. Finally, only data collected during the period from 24 to 38 weeks were taken into analysis due to the lack of AOT observation during week 39, so the analyzed data covered the period from 13 June to 18 September 2016.

Diurnal Patterns of VPD
The NEP vs. SWin dependency was used during this analysis and since carbon dioxide uptake of the peatland has a regular diurnal pattern, the hours of day values (UTC + 1) were used as a proxy for solar angle [100] (Figure 4 and Figure 5). At the Rzecin site, the reduction of NEP was found between 12PM and 15PM ( Figure 4) and it is associated with VPD values higher than 20 hPa. This observation corresponds with literature findings where VPD threshold values between 15 and 20 hPa were reported [96,101,102].   Thus, the data collected during the periods when VPD exceeded the threshold of 20 hPa was excluded from further analysis. The result of this selection has been presented in Figure 5d.

Diurnal Patterns of Micrometeorological Parameters
The meteorological variables of Ta, VPD, DI, and NEP were presented in the context of both local time (solar angle proxy) and SWin ( Figure 5). The afternoon values of Ta and VPD are higher than those observed in the morning (Figure 5a,b). The diffusion index (DI) is also higher in the afternoon than in the morning (Figure 5c). There is a drop of DI values between the hours of 15PM and 18PM, but the values are still higher than in the analog time of the morning and it corresponds to slight Ta increase that will be determined by higher insolation during this part of the day (Figure 5a). The ecosystem was more effective in terms of carbon dioxide absorption (higher NEP) during the afternoon than in the morning hours (Figure 5d).

ART Model Parametrization
In order to estimate the direct and diffuse solar irradiance using photometric measurements, the ART model was used. The estimated values of solar irradiance (SWmod) were compared with the measured ones (SWin). There were very good matches found for total and diffuse irradiance, where b1 = 1.06, b0 = 60.3, R 2 = 0.913 (Figure 6  The ART overestimates some of the shortwave flux density values (SWmod) ( Figure 6). This is probably due to the fact that the instantaneous values of AOT are applied for the 30minute shortwave radiation flux estimation while pyranometer's observations are carried out over the whole observation period. Additionally, both the field of view and measuring procedures of those sensors differ substantially and it can lead to SWmod overestimation during the moments when sparse clouds are present in the atmosphere.

EP Model Parametrization
The EP model was parameterized on the basis of measured values of Ta, total (I t ), direct (I r ), and diffused (I f ) irradiance. Nonlinear least-squares model was applied and all model parameters were found to be statistically significant (MAE = 2.36 µmol·m −2 ·s, RMSE = 3.09 µmol·m −2 ·s, NRMSE = 8.7%, N =1 034 (Table 1).

Effects of Meteorological Conditions and Optical Properties on NEP
Three levels of DI values were used during these calculations, where the observed DI values were increased and decreased by 0.1 (Table 2, Figure 7). This range of DI values is considered typical for cloud free conditions [47].  The values of the measured DI that were applied in the calculation ranged between 0.16 and 0.57 during the analyzed day. These values are typical for sunny conditions, however, the additional turbidity virtually introduced by the modification of DI can be determined by the presence of clouds, e.g., thin cirrostratus type or increased transparency of air. Thus, only the simulations performed with the clear sky ART model enables the estimation of the cloud free atmosphere optical parameters impact on the peatlands production.
The conditions from the 23rd of June 2016 were applied again as input data in the EP model, however the ART model was used for the calculation of the total, direct, and diffused PPFD. These estimations were done for three AOT500 values: 0.09, 0.13, and 0.17 that are the 1st quartile, median, and 3rd quartile of AOT500 population, respectively (Table 3, Figure 8). The change in diffusion index (DI) caused by the addition or subtraction of 0.1 from the measured values induced an NEP change of 0.68-0.98 µmol·m −2 ·s −1 (6.1-42.3%) and −1.55 to −0.79 µmol·m −2 ·s −1 (−49% to −10.5%) of CO 2 uptake, respectively ( Table 2). The application of the ART model allowed for the assessment of aerosol impact on peatland NEP. The increase of AOT500 value from 0.09 to 0.13 and 0.17 determined the rise of NEP in the range between 0.08 to 0.22 µmol·m −2 ·s −1 (1.9-3.1%) and 0.16 to 0.48 µmol·m −2 ·s −1 (3.4-5.7%) of NEP, respectively (Table 3). On the basis of these results, it was found that the method is responsive to AOT changes and there is a potential in CO 2 uptake that is noticeable in the AOT increase from 0.09 to 0.17.
Aerosols can both absorb and scatter the radiation, thus the increase of AOT attenuates the total irradiance (I t ) and increases the share of diffuse radiation (DI) in this energy flux. Usually, the increase of DI determined by AOT500 increase results in ecosystem production (NEP) rise, however during the morning period the photosynthesis reduction determined by I t attenuation prevails the gain obtained due to diffuse radiation impact. This situation is observed under low I t values (<250 umol·m −2 ·s −1 ) when NEP values are close to zero (gross ecosystem production compensates ecosystem respiration) (Figure 8).

Discussion
Peatlands, due to their hydrogenic origin, are very sensitive to shifts of both thermal and hydrological conditions, but seem to be sensitive also to diffuse radiation that can determine their productivity rise as a result of higher light use efficiency, which is known as a diffuse fertilization effect [103]. The application of locally parametrized ART and EP models allowed the estimation of the peatland's productivity change as a result of the radiation scattering modifications determined by AOT change. Such a strategy of study does not often appear in the context of the research of diffuse radiation impact on terrestrial vegetation productivity [104].
The short-term afternoon reduction of NEP (Figure 4) was observed in Rzecin. This limitation of NEP is probably determined by both the high leaf/air temperature that determines the heat stress, as well as the high vapor pressure deficit that reduces the gas exchange between leaf and ambient atmosphere [46]. This midday depression of the net ecosystem photosynthesis can be reduced by the aerosol presence in the atmosphere due to the leaf temperature decrease, which was also found in studies of terrestrial vegetation in East Asia [105].
The initial analysis of the collected data showed that NEP is noticeably reduced during periods when high values of VPD are found. Similar strong negative sensitivity to VPD of the Siberian taiga NEP was found on the basis of cloud and fire events analysis [106]. Thus, the stepwise data elimination showed that the afternoon NEP depression was not observed when VPD was lower than 20 hPa. This value corresponds with the results obtained by other authors for spruce, semi-arid steppe, and alpine shrub vegetation [16,107]. The wetland plants (well-watered environment) react to VPD similarly to other non-aquatic plants, and this indicates that the stomata of all plants act according to a common scheme.
The diurnal patterns of Ta and VPD have characteristic shapes where the afternoon values of both have been found to be higher than in the morning (analogous hours). Maximum values of both Ta and VPD were found around noon. The temperature rise induces CO 2 efflux (ER). Additionally, higher air temperature results in an increased VPD that may inhibit photosynthesis (GEP). Thus, the higher values of Ta and VPD observed during the afternoon seemingly negatively affect NEP. In contrast to those two parameters, afternoon NEP values are higher than those observed during analogous hours in the morning. Response of DI between 15PM and 18PM can be unexpected (Figure 5c). This afternoon's DI decrease is probably the effect of the convection process disappearance that leads to a situation where the sun's disc is less often obscured by the clouds during low solar angles. However, in general, the DI values in the afternoon are still higher than in the morning, and this parameter seems to be a factor that is able to compensate and even surpass the reduction of NEP determined by higher Ta and VPD. The results of this analysis are in very good agreement with other authors findings in case of Beech forest [100].
The Canadian peatland study showed no evidence of higher net ecosystem CO 2 exchange under diffuse light conditions; however, this study area was covered mainly by the vegetation characterized by noticeably lower LAI (maximum 1.3 m 2 ·m −2 ) [50].
Typical values of DI for cloud free conditions are usually found in the range between 0.1 and 0.3 [16]. In this analysis, real measured values (0.16-0.57) were used as an input for NEP estimation using the EP model. The impact of a diffusion intensity change was estimated by the application of both the increase and decrease of DI by 0.1. Those modifications resulted in the change of DI values and its minimum values remained in a reasonable range from 0.06 to 0.26. For all I t values, the NEP rose by 0.68-0.98 µmol·m −2 ·s −1 (6.1-42.3%) for a DI increase of 0.1 and declined by 0.79-1.55 µmol·m −2 ·s −1 (−49.0% to −10.5%) for a DI decrease of 0.1 ( Table 2). The decrease of DI by 0.1 results in higher NEP reduction (−49.0 to −10.2%) than a DI increase (6.1-42.3%). Such asymmetric reactions of the ecosystem to changes of DI values can be explained by the lower radiation use efficiency of upper leaves that gain the radiation load under less scattered radiation. The same EP model was used for the estimation of NEP during the cloud free conditions when the AOT500 value increased from 0.09 to 0.13 and from 0.13 to 0.17 ( Table 3). The direct effect of AOT increase resulted in a reduction in incoming radiation, but simultaneously it determined an increase of NEP for I t values higher than 300 µmol·m −2 ·s −1 . The estimated NEP rises were determined by the changes in AOT500 value from 0.09 to 0.13 and from 0.13 to 0.17 resulting in an NEP increase in the range of 0.08-0.26 µmol·m −2 ·s −1 (1.9-3.1%) and 0.15-0.48 µmol·m −2 ·s −1 (1.6-2.7%), respectively (Table 3). In the case of I t values below 300 µmol·m −2 ·s −1 , a reduction of NEP was found, and under these conditions the increase of diffuse radiation did not compensate for the loss of NEP caused by I t decrease. The presented results correspond with other studies (e.g., [108][109][110]) where the effects of atmospheric aerosol presence on the absorption of CO 2 by terrestrial ecosystems showed opposite tendencies. An upward trend under moderate aerosol concentrations and cloudless conditions was found when the Mount Pinatubo eruption enhanced the deciduous forests photosynthesis rate by 23% and 8% in 1992 and 1993, respectively [44], while the increase in aerosol loading by 0.09 to 0.16 caused the boreal and hemiboreal forest GPP to rise by 6-14% [111]. The opposite situation was found during the significant reduction of total radiation determined by high AOT values, i.e., higher than 2.7 [110], and the forest canopy modeling study showed that crossing the AOT threshold of 1.0 resulted in canopy photosynthesis rate reduction despite the continuous light use efficiency increase [112].
These results are not surprising, since the diffuse irradiance values obtained from the ART model (cloud free conditions) were at least 40% lower than the observed one (b1 = 0.5948, b0 = 26.1594, R 2 = 0.667). The agreement between simulated and measured diffuse irradiance values was also reduced (lower R 2 ) by the possible presence of clouds on the sky dome ( Figure 6). In other words, the disagreement between simulated and observed values may be rooted in uncertainty of the cloud free conditions described in the results section.
The changes in NEP determined by both DI and AOT500 modifications are weaker for higher values of I t . This result indicates that ecosystem production potential for higher I t is usually limited due to reaching the maximum capacity of plants to capture light and fix CO 2 [104]. However, even in these circumstances the activation of lower leaves can determine higher CO 2 uptake of the canopy.
These findings show that NEP of short plant canopies is sensitive to diffuse radiation. This indicates that a complex spatial structure of a vegetation layer is the primary feature, more than plant size, causing the interaction of ecosystems with diffused radiation. The results presented here show that aerosols can noticeably impact the peatland CO 2 uptake. This determines the necessity of peatland-atmosphere interaction studies at larger spatial and temporal scales to account for aerosol variability. This extrapolation in time and space can be obtained by the application of observational networks that operate at a global scale. The following networks, then, can be the source of data for such extrapolation: integrated carbon observation system (ICOS) [113], baseline surface radiation network (BSRN) [114], and aerosol robotic network (AERONET), [37] where the ecosystem production, solar irradiance, and aerosol optical thickness are permanently measured, respectively.
Changes in the diffuse solar radiation (SWdif) received at the Earth's surface have been observed over the last several decades [18,[20][21][22]115]. In detail, changes in SWdif vary over the world, depending on the geographical location and measurement period. It has been noticed that there was a decrease in the annual mean series of SWdif in Spain from 1985-2010 of 2.1 W·m −2 per decade [116]. These results agree with the trends reported since the 1980s from a rural site in Estonia [19]. However, in the period from 1955-1992, no significant trend of SWdif was found there. A weak tendency towards a decrease in diffuse solar radiation during the period 1955-2005 was obtained in Moscow [20]. The decline in SWdif by 2.44 W·m −2 has also been observed in Germany from the 1960s to the late 1990s. The most likely causes of increases in global irradiance are the long-term decrease in the number of aerosols and longer sunshine duration [116].
The observed decrease in SWdif is in line with the brightening period reported in many regions of the world, which implies less scattering due to the decrease of aerosols and/or cloud cover during this period [117]. In the observed trends in SWdif, there are some exceptions. For example, Stanhill [18] found a statistically significant increase in SWdif at Dublin airport station. There are also some places where positive trends of SWdif were observed, e.g., U.S. (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) [118] and India   [22].
The diffuse solar radiation study in Poland shows the spatial and temporal variability of trends, and from the early 1970s to the mid -1990s, the annual DI mean series showed a significant increase in Warsaw [119], Wrocław, the Upper Silesian Industrial Region, and the Rybnik Coal Area [120], while a decrease in SWdif was found at the beginning of the 1980s at Wrocław station [121]. Generally, the decrease of the scattering process intensity prevails and the reduction of the ecosystem production due to DI decrease should be observed on a global scale [14]. However, the reaction of peatlands still requires estimation at the global scale. There is also the question of the reaction of peatland net production to a substantial increase in diffuse radiation determined by stratospheric eruptions, e.g., Mount Pinatubo [122].
The applied ART model estimated total solar radiation in a very realistic way, while the estimation of diffuse radiation flux is underestimated. This was mainly determined by the assumption of cloud free conditions. The measured diffuse radiation was found to have higher values than the modeled one and it was the effect of cloud presence. Thus, the application of the ART model gave us the opportunity to estimate the diffuse radiation impact on the cloud free moments and allowed for a separate study of aerosol impact on the peatland CO 2 uptake. EP model was parameterized for conditions when VPD values were lower than 20 hPa and it is a limitation of the presented methodology. However, this data selection was necessary to omit the most important non-radiative limiting factor of plant photosynthesis. However, the obtained results show the the aerosol observations [35] in the atmosphere (their lifetime, properties and long-range transport) are important for the peatland or other ecosystems' CO 2 uptake estimations.
The combination of the EP and ART models enabled the impact of aerosols on Rzecin peatland CO 2 uptake to be deconvolved. This approach provides the linkage between the ecosystem production and optical properties of the atmosphere, and allows for the shift of the scale of the study by application of remotely sensed data. The impact of diffuse radiation is usually considered an important factor that determines the NEP of higher plant canopies. This has been the primary reason that research related to diffuse radiation impact on ecosystem production has focused on forest ecosystems [14,16,[42][43][44][45][46][47][48]. This study suggests that shorter vegetation investigations should be carried out more intensively, as fewer studies have been conducted on shorter plant canopies [43,49,50]. The presented analyses enrich our knowledge about the impact of aerosol presence on central Europe peatland CO 2 exchange, which provides the basis for further and deeper considerations in this context. This analysis is also interesting due to the future cloud behavior simulations. It is predicted that the occurrence of clouds in some regions of the world will decrease [123] or even that some of their types (such as marine stratocumulus) will disappear [124] along with an increase in the concentration of greenhouse gases in the atmosphere. Moreover, even if the impact of aerosols on the peatland carbon balance is transient, further work will be required to characterize and quantify the diffuse radiation impact on the CO 2 exchange of these valuable ecosystems.

Conclusions
The complex character of the interaction between the ecosystems and atmospheric conditions required the application of models that provided the basis for the deconvolving of single factors, e.g., aerosol and/or cloud impact on ecosystem production. The study focused only on the impact of optical parameter changes, determined by the aerosol presence in the atmosphere, on peatland net CO 2 uptake during cloudless conditions.
Presented research elucidates the interactions between the optical parameters of the atmosphere, (e.g., DI and AOT500) and peatland production. Moreover, it clearly indicates that the scattering radiation process determined by the presence of aerosol modifies the amount of absorbed CO 2 from the air.
These findings indicate a non-negligible impact of the optical properties of the atmosphere on the peatland NEP, thus the scattering process should be taken into account during further mass and exchange study of sphagnum dominated peatlands. It also indicates that future changes of diffuse irradiance will modify the CO 2 uptake of transitional peatland. The study reported here is ont a 'local' scale. However, it indicates that more work needs to be done regionally and globally to better understand the influence of direct and diffuse solar radiation, particularly on peatland ecosystems and the accumulation or emission of carbon for climate change studies, if these studies are to be robust and included in climate change models.