The trend of the oxidants in boreal forest over 2007-2018: comprehensive modelling study with long-term measurements at SMEAR II, Finland

Major atmospheric oxidants (OH, O3 and NO3) dominate the atmospheric oxidation capacity, while H2SO4 is considered as a main driver for new particle formation events. Although numerous studies have investigated the long-term trend of ozone in Europe, the trend of OH, NO3 and H2SO4 at specifc sites are to a large extent unknown. In this study, we investigated how the trends in major atmospheric oxidants (OH, O3 and NO3) and H2SO4 changed in southern Finland during the past 12 years and discuss how these trends relate to decreasing emissions of regulated air pollutants in Europe. 1


INTRODUCTION
Understanding the atmospheric oxidants (OH,O 3 and NO 3 ), their reactions and related processes is important as they are the main "cleaning protagonists" of the atmosphere. Many trace gases, such as methane (CH 4 ), volatile organic compounds (VOCs), nitrogen oxides (NO x = NO + NO 2 ) and sulphur dioxide (SO 2 ) are removed from the atmosphere by oxidation reactions. During the day, the hydroxyl radical (OH) is the dominant oxidant produced by photochemical processes in the troposphere al., 2008;FAO Global Forest Resources Assessment, 2015) and boreal vegetation is dominated by evergreen coniferous trees that produce signifcant amounts of biogenic VOCs (BVOCs), mainly isoprene (C 5 H 8 ), monoterpenes (C 10 H 16 ) and sesquiterpenes (C 15 H 24 ) (Hakola et al., 1998(Hakola et al., , 2006Rinne et al., 2009). Studies of the OH-reactivity in forest canopies have suggested large emissions of unknown reactive BVOCs (Mogensen et al., 2015;Praplan et al., 2019).
Studies on long-term trends of oxidants can provide an insight on how the atmospheric oxidation capacity evolves against the background of climate and local changes in the environment. Several studies have investigated the trends of atmospheric oxidants in Europe. In their studies, Wilson et al. (2012) and Yan et al. (2018) showed a general decreasing trend in ozone concentrations due to the decrease in NO x -emissions. Numerous studies have investigated global OH trends using chemical transport models or retrieval of remote sensing of methylchloroform (CH 3 CCl 3 , MCF) (e.g. Montzka et al., 2000;Prinn et al., 2001, Kirschke et al., 2013. Montzka et al. (2011) found a small interannual OH variability, indicating that global OH is generally well buffered against perturbations. In situ long inter-annual OH measurements are relatively rare. However, one study by Rohrer et al. (2006) showed that there was no trend in the level of OH in the Hohenpeissenberg data set during the studied period 1999-2003 (estimate annual trend to be less than ± 2.5% yr -1 ) and that there was a positive correlation (r=0.941) between OH and the photolysis frequency of ozone, J(O 1 D). Long term trends of NO 3 are rarely studied, and only a few modelling studies on the long-term NO 3 trend exist (e.g. Heintz et al., 1996).
A considerable number of feld campaigns, in which OH concentrations were measured, have been compared to the results of modelling simulations (e.g. Eisele and Tanner, 1991;Holland et al., 1995;Petäjä et al., 2009). Most modelling studies reproduced the OH concentration within the uncertainty range of the OH measurements, including clean (e.g. Tan et al., 2001;Ren et al., 2005;Kubistin et al., 2010;Dlugi et al., 2010;Kanaya et al., 2012;Regelin et al., 2013) and urban areas (e.g. Heard et al., the literature to overcome this problem like the "Uniform Fill-In" or the "Log Fill-In" methods discussed and tested by Cohen and Ryan (1998). However, as all data below LOD are unknown no method predicts their distribution correctly which makes it difficult it difficult to choose a single technique that will be best at all times for various parameters. In Table S1 in the supplementary material we calculated the amount of data points for SO 2 and NO above LOD (the two parameters with the highest amount of data below LOD) for different percentile ranges in each year to investigate if a trend in the below LOD data exists.
The measured CH 4 concentrations in 2014 were used as input in SOSAA for the year 2014. For other years, an annual growth rate of 6 ppb yr -1 was assumed, and the input CH 4 concentrations from 2014 were thus added or subtracted a multiple of the annual growth rate depending on the year difference.
The growth rate were chosen from the 'NASA Earth Observatory' website and represent the methane increase for 2007-2013 (https://earthobservatory.nasa.gov/images/87681/a-global-view-of-methane).
The condensation sinks (CS) for H 2 SO 4 and HNO 3 were provided as an input for the model. The CS was calculated based on the particle size distribution measured by a DMPS (particle diameters 3-1000 nm) and an APS (particles with aerodynamic diameters 0.5-20µm) system (Pirjola et al., 1998;Kulmala et al., 2001), and the hygroscopic growth effect was corrected based on Laakso et al. (2001).
Similarly to the meteorological input data, the input mixing ratios and the CS were also linearly interpolated to 60 s time resolution to match the simulation time step of the emission and chemistry modules.

STATISTICAL METHODS
The daily/daytime/nighttime trends of variables were calculated based on their daily/daytime/nighttime mean or median values. Whether to use mean or median for a variable is determined by its data value distribution. If the data are logarithmically distributed (O 3 , CO, CS, EM-MON, MON, OH, HO 2 , H 2 SO 4 , NO 2 , N 2 O 5 , and NO 3 ), the median values are used. Here we should notice that although the data value distributions of SO 2 and NO are also logarithmically distributed, we still used their mean values.
The reason is that more than 50% of their measured concentrations lie below the LOD, which results in that their median values are equal to LOD. For other variables (temperature, RH, and solar irradiance), the mean values were used. For the logarithmically distributed variables (besides the variables mentioned above, SO 2 and NO are also included here), the daily/daytime/nighttime linear trend fttings were conducted on the logarithm with base 10 of their respective median or mean values. For other variables (temperature, RH, and solar irradiance), the linear trend fttings were performed directly on their respective mean values.
Bootstrapping was used to estimate the confdence interval of the trend (Wilks, D., 1997, Asmi at al., 2013. We frst ftted a linear trend to the time series and created a new data set by taking random samples from the original residuals (differences of the data values and the ftted linear trend) and adding these to the linear part. Then a new linear ft was made to this new data set. This procedure was repeated several times (typically 1000 iterations). Here the idea was to test the monotonicity of the trend. The smaller the differences in the ftted trends were after many such iterations, the original trend was more likely to be monotonic. To get the confdence interval we examined the 5th to 95th percentile range of the slopes obtained from bootstrapping iterations: if all of the slopes in this range were either positive or negative (thus not containing zero trend), we concluded that the likelihood of there being a trend was higher than 95% (p<0.05) and thus statistically signifcant. To get another estimate of the monotonicity of the trends we also used the Mann-Kendall test for autocorrelated seasonal data (Hipel and McLeod, 1994;Hussain et al., 2019), and the p-values are reported in Table 1 under P MK . The MK test is more conservative, but both our tests agree in the sense that wider confdence interval or larger p-value indicate larger yearly variation and hence the prognostic capacity of the trend is smaller. Finally, the relative changes (and the 90% confdence interval from the bootstrapping test) which are shown in Table 1 are in linear scale for all variables, describing the actual change in % yr -1 or variable units yr -1 . The average trendumber feld. You see the Edit Fields dialog. shown in time series plots is obtained with a 1-year running median (window size of 182 days), see Ma et al.(2016) for a detailed description.

RESULTS
The results will be presented in 6 subsections: (1) a short validation of the meteorological module, (2) the trends of measured gases, (3) BVOCs (observation and model inter-comparison), (4) trends and campaign model-observation inter-comparisons for the main oxidants (O 3 , OH and NO 3 ), (5) trends and campaign model-observation inter-comparisons for sulphuric acid and (6) comparisons of proxies for OH, NO 3 and H 2 SO 4 with the model results.  Here the daytime is defned as the time period between one hour after sunrise and one hour before sunset; nighttime is defned between one hour after sunset and one hour before sunrise and the daily values are averaged for 24 hours. In the following subsections we will discuss the values for single parameters from Table 1 in more detail.

METEOROLOGICAL DATA ANALYSES
Meteorology is one of the major drivers for the change in atmospheric composition. We compared several measured meteorological parameters with the model outcome to validate the performance of the meteorological module in SOSAA. While temperature, water vapour and wind speed were nudged with the measurements, the heat fluxes and net radiation were simulated and their comparison with measurements provide an insight into the simulated energy balance above the forest canopy.   Here, the net radiation is calculated as the total incoming short-and long-wave radiation subtracting the total outgoing short and long-wave radiation at the canopy top.
Thus positive values represent more incoming than outgoing radiation and vice versa. The modelled daytime values agree well with the measurements in autumn (September, October and November) and early winter (December and January), while the model overestimates the measurements from late winter (February) through summer, with the exception of July. The overestimation occurs mostly at noon with the averaged noon peak overestimated values ranging from ~20 Wm -2 to ~100 Wm -2 . By contrast, the modelled nighttime net radiation underestimated the measurements by about ~10 Wm -2 tõ 50 Wm -2 on average from September to December and from January to March. In general, the model is consistent with the measurements, and is able to capture the diurnal pattern and seasonal trend of net radiation above the canopy. Therefore, considering the simulation results of SH and LH discussed above, the model can predict a reasonable energy balance inside and above the canopy. In Table 1, the temperature and relative humidity represent the analyses of the measured data which are used to nudge the model as mentioned above. The trend for daily mean temperature shows an increase

TREND OF INORGANIC GASES AND CS
The main inorganic gases (CO, O 3 , NO, NO 2 and SO 2 ) that are read in as input to SOSAA reflect the influence of human impact on a regional scale. Carbon monoxide for example has a lifetime of approximately 1-3 months (Seinfeld and Pandis, 2006), and the concentration levels reveal the impact of large regional to hemispherical features. Nitrogen oxides and sulphur dioxide have lifetimes of days to weeks, respectively, and they are mainly related to local or regional changes. At a clean background station like SMEAR II, their concentrations are often below the LOD of the instruments.
The 12-years concentrations of fve measured trace gases (CO, O 3 , NO, NO 2 and SO 2 ) and the aerosol condensation sink (see Table 1 and Fig. S2

LONG TERM TIME SERIES OF MONOTERPENES
Previous studies used the empirical proxies method to investigate monoterpene seasonal and diurnal variations (Kontkanen et al., 2016), which may contain high uncertainties. Based on the long-term and evaluated simulations of monoterpenes concentrations by SOSAA, we analysed the long-term trend of monoterpenes concentrations at SMEAR II and presented the results in Ozone is one of the most important oxidants in the atmosphere and thereby its trend during 2007-2018 will be discussed here based on the continuous measurements at SMEAR II. Fig. S2b and Table 1 show that the ozone concentrations are relative stable over the 12 years (change is -0.11 (

HYDROXYL RADICAL -OH
The hydroxyl radical, OH is the most important oxidant in the troposphere and it is the major "cleaning protagonist" in the atmosphere by reacting with nearly all trace gases including the vast number of VOCs emitted from the boreal forest. OH is also the most important sink term for methane (CH 4 ), the second most important greenhouse gas, responsible for approximately 20 % of induced global radiative forcing since pre-industrial times (Turner et al., 2018). OH is also crucial for the sulphuric acid production and in this way, it indirectly influences the formation of secondary organic aerosols (see subsection 3.5). Therefore, it is important to study the trend of the OH concentrations and to investigate whether increased temperature or changes in the gas-phase composition in the atmosphere during the last decade at the SMEAR II station had an impact on the OH concentrations.
The OH measurements are difficult and expensive and therefore measurements of OH at the SMEAR II are rare. In this study, the measurements from the EUCAARI 2007  and the COPECC-HUMPPA in 2010 (Williams et al., 2011) campaigns were used to evaluate the model performance. To test the simulated OH concentration, we compared measured data from these two campaigns against the model results (Fig. 5).   The increasing trend of the OH concentrations at the SMEAR II station is somehow surprising on the frst view, considering the increase of the monoterpene concentrations for the same period. Besides reacting with OH and being a sink term for OH, monoterpenes also produce OH through reaction with ozone. This is the main source of OH during dark conditions at SMEAR II and the nighttime increase of the hydroxyl radical by about 2.8% yr -1 is related to this mechanism. However, the absolute OH nighttime concentrations are less than one tenth of the daytime values and therefore the ozonolysis of monoterpenes has only a small contribution to the daytime increase of OH in our calculations. As pointed out in earlier publications (e.g. Boy et al., 2006;Praplan et al., 2019), carbon monoxide (CO) is the main sink term for OH and accounts for about 40% of the removal of OH in the troposphere. CO has a lifetime of 1-3 months (Seinfeld and Pandis, 2006) and has decreased since the 1990s, as shown by Greenland frn air records (Wang et al., 2012;Petrenko et al., 2013) and by surface flask samples collected at few sites (Khalil and Rasmussen, 1994;Novelli et al., 2003;Gratz et al., 2015;Schultz et al., 2015). The drop of CO in our study (0.5% yr -1 ) affects the OH trend more than the increase of BVOCs as other important parameters concerning the OH production (e.g. solar irradiance or ozone concentrations) are more or less unchanged during the 12 years. Additionally, the hydroperoxyl radical (HO 2 ) surprisingly indicates a positive trend of +2.89 (+0.81; +4.74) % yr -1 , even stronger compared to OH. This is related to the decrease of the NO concentrations as the nitrogen monoxide is by far the most important sink term for HO 2 (Boy et al., 2006).

NITRATE RADICAL -NO 3
The nitrate radical simulated in this study can't be validated by observations, due to lack of measurements at the SMEAR II. However, by constraining SOSAA with accurately measured NO 2 concentrations, we assume that the predicted NO 3 concentrations are reasonable. The rapid photolysis of NO 3 and the reaction with NO typically reduces its lifetime to a few minutes during daytime. The main contribution of the nitrate radical to the oxidation capacity of the atmosphere is during nighttime.
Based on this, we will focus our analysis on the nighttime period.  The reason for this pattern which is visible in the model outcomes for all years is a combination of mainly three effects: SO 2 , one of the two main precursors for H 2 SO 4 , is peaking in late winter and early spring and OH reaches its yearly maxima in spring. Additionally, the condensation sink, representing the rate of how fast sulphuric acid molecules will condense on the existing particles, has a clear maxima in summer (SF2 in supplementary material). These three parameters are mainly responsible for the sulphuric acid pattern. Note that a similar pattern has been observed for the occurrence of NPF events at the SMEAR II for several years (Nieminen et al., 2014). Nieminen at al. (2014) predicted the trend of sulphuric acid based on a proxy calculation (see next subsection) with -1.3% on NPF days and -0.3% on non-NPF days per year for the years 1997-2012. In our study, we applied SOSAA simulations for the years 2007-2018 for the same location. Our model results predict a stronger decrease of daytime H 2 SO 4 of -5.12 (-11.39; -0.52) % yr -1 (see Table 1).
However, the confdence interval of this trend is quite broad, which tells that caution should be taken to interpret this trend too far to the future. The trend in the studied time span is greatly influenced by the large yearly variation, especially year 2016 deviates, which can be seen in Fig. 9 During the last years, several proxies have been developed for compounds like the hydroxyl or the nitrate radical, due to absence or sparse long-time observations for these parameters. In this subsection, we will compare some of these proxies with the outcome of our model simulations for SMEAR II. This, however, should not be seen as a validation of the proxy but rather to investigate how well simulations agree with them. The proxies compared were developed based on datasets from the SMEAR II. The frst proxy we compare is for the OH radical. It is based on a publication by Nieminen at al. (2014). The results of the proxy, together with the outcome of the model simulations, are presented in Fig. 10. In this proxy, the hydroxyl radical is calculated as; [OH] = ((8.4 * 10 -7 / 8.6 * 10 -10 ) * UVB 0.32 ) 1.92 Eq. 1 Here UVB presents the ultraviolet irradiance measured at SMEAR II. The modelled and the proxy OH concentrations show a similar trend from October to April but start to diverge from May to September with the highest discrepancy around late May to late June. During this time the modelled OH concentrations are about two-fold higher compared to the proxy. Later in the year from July to September the modelled values are still higher compared to the proxy but the discrepancy decreases.
As pointed out in subsection 3.4.2, we believe, based on OH-reactivity measurements at SMEAR II, there exist missing compounds reacting with OH. Taking this into account and assuming that the missing compounds originate from the local ecosystem with maximum emissions during the most biologically active period, the modelled OH concentrations are potentially too high during spring and summer and the proxy would be more accurate at these times. However, as long as these unknown compounds are not identifed and no long-term measurements of OH at the SMEAR II exist, any fnal conclusion about whether the proxy or the model is more correct can only be speculation. The next proxy we compare with our model simulations is for sulphuric acid. It is based on a parameterization method from Petäjä et al. (2009). The proxy is calculated as follows,

[H2SO4] = k * [SO2] * UVB / CS Eq. 2
Here UVB stands for the ultraviolet irradiance, CS for the condensational sink and [SO2] for the gas phase concentration of sulphur dioxide. The scaling factor k is an empirically derived factor, which scales the proxy variables to correspond to the measured sulphuric acid concentrations. As already pointed out in the previous subsection, the modelled sulphuric acid concentrations show a clear peak in early spring and then a nearly continuous decrease for the rest of the year (see Fig. 11). The proxy follows this pattern almost identically from April to September but exceeds the modelled data approximately by a factor of 2 during the months from October to April. As the measured values for this period seems to agree well with the model results, we conclude that the proxy in autumn and winter overestimates the H 2 SO 4 concentrations considerably and we assume that for this parameter, the model provides a more realistic picture of the sulphuric acid concentrations.  The last proxy we want to compare against model results is for the nitrate radical. It is based on a publication by Kontkanen et al. (2016). The concentration of NO 3 is calculated based on the following equation: Here kO3+NO2 is the temperature-dependent reaction rate coefficient between NO 2 and O 3 , which was calculated from a temperature-dependent relation (Atkinson et al., 2004; see Table A1 in Kontkanen et al., 2016). τ NO3 is the lifetime of NO 3 and a detailed description on the prediction of τ NO3 is available in the manuscript by Kontkanen et al. (2016).
Besides O 3 , the nitrate radical is the most important oxidant during nighttime and has an important contribution to BVOC oxidation during nighttime (Mogensen et al., 2015). However, until now, NO 3 measurements at SMEAR II are rather limited and most of the existing data achieved during the COPECC-HUMPPA campaign in 2010 (Williams et al., 2011) and the IBAIRN campaign in 2016 (Liebmann et al., 2018) were below the LOD.  Thus, a robust estimate for NO 3 concentration for the whole year is needed. Previous studies assumed a steady state between the production of NO 3 from the reaction between O 3 and NO 2 and the removal of NO 3 (Kontkanen et al., 2016). This gap could be flled by our study. By applying this approximation, we derived the monthly median NO 3 proxy between 2007-2018. A comparison between the proxy and the model data shows that the long term trend from both methods are in very good agreement (Fig.   12). As those two methods applied here agree very well, it is likely that the predicted values for NO 3 are reliable and could be applied in further studies.

SUMMARY AND PERSPECTIVES
In this study we investigated the trends of various measured and modelled meteorological parameters   with the OH. Therefore, the predicted OH trend shows that the climatic temperature increase (~0.8K in 12 years at SMEAR II) and the followed rise in BVOC emissions is buffered by a decline of carbon monoxide. In case the current negative trend in CO continues -mostly related to improved combustion techniques -the OH will slightly rise or at least stagnate at the present level which reflects a positive impact on the atmospheric oxidation capacity. Vice versa is the situation for NO 3 showing a nighttime decrease by -3.92 (-6.49; -1.79) % yr -1 which is caused by the drop of nitrogen dioxide. As all anthropogenic NO x emissions in Europe have decreased signifcantly during the last decades (EAA report No 12/2018) and are predicted to decrease further, we expect that the nitrate radical will continue to drop in the future.
Sulphuric acid was investigated as it is one of the most important precursors of new particle formation (NPF). The outcome of our study indicate that the sulphuric acid concentration is decreasing with -5.12 (-11.39; -+0.52) % per year during daytime, which likely is related to the reduction in the emissions of sulphur dioxide in Europe during the last decades (EAA report No 12/2018). In case the negative trend of sulphuric acid (steered by SO 2 ) will continue in the next decades, it could affect the amount of NPF events in the boreal region signifcantly. However, whether or not this will have a positive or negative impact on our future climate is currently unclear. In the past, it was typically assumed that NPF events will provide more CCN followed by more cloud droplets, leading to an increased albedo through "brighter" clouds (e.g. Makkonen et al., 2012). In this way, NPF will cool the planet and counteract the effect of greenhouse gases. However, quite recently Roldin and co- sulphuric acid should be seen as positive or negative. However, it is certain that H 2 SO 4 has decreased in the last decades and, most likely, will continue to drop in the future.
Proxies are commonly applied in case limited amount of parameters are measured and no detailed model simulation are available. We compared concentrations for OH, NO 3 and H 2 SO 4 calculated from proxies (Nieminen at al., 2014;Petäjä et al., 2009;Kontkanen et al., 2016) with our model outcomes.
Our comparisons showed that the proxies for the OH and H 2 SO 4 agree at certain times of the year very well with the model results but also differ signifcantly during other periods. For the nitrate radical, the model and proxy results are in good agreement. Competing interests. The authors declare that they have no conflict of interest.

Special issue statement.
This article is part of the special issue "Pan-Eurasian Experiment (PEEX)".
It is not associated with a conference.