Future vegetation – climate interactions in Eastern Siberia : an assessment of the competing effects of CO 2 and secondary organic aerosols

Disproportional warming in the northern high latitudes and large carbon stocks in boreal and (sub)arctic ecosystems have raised concerns as to whether substantial positive climate feedbacks from biogeochemical process responses should be expected. Such feedbacks occur when increasing temperatures lead, for example, to a net release of CO2 or CH4. However, temperature-enhanced emissions of biogenic volatile organic compounds (BVOCs) have been shown to contribute to the growth of secondary organic aerosol (SOA), which is known to have a negative radiative climate effect. Combining measurements in Eastern Siberia with model-based estimates of vegetation and permafrost dynamics, BVOC emissions, and aerosol growth, we assess here possible future changes in ecosystem CO2 balance and BVOC–SOA interactions and discuss these changes in terms of possible climate effects. Globally, the effects of changes in Siberian ecosystem CO2 balance and SOA formation are small, but when concentrating on Siberia and the Northern Hemisphere the negative forcing from changed aerosol direct and indirect effects become notable – even though the associated temperature response would not necessarily follow a similar spatial pattern. While our analysis does not include other important processes that are of relevance for the climate system, the CO2 and BVOC–SOA interplay serves as an example for the complexity of the interactions between emissions and vegetation dynamics that underlie individual terrestrial processes and highlights the importance of addressing ecosystem–climate feedbacks in consistent, process-based model frameworks.


Introduction
Warming effects on ecosystem carbon cycling in northern ecosystems (Serreze et al., 2000;Tarnocai et al., 2009) and the potential for large climate feedbacks from losses of CO 2 or CH 4 from these carbon-dense systems have been widely discussed (Khvorostyanov et al., 2008;Schuur et al., 2009;Arneth et al., 2010).Other biogeochemical processes can also lead to feedbacks, in particular through emissions of biogenic volatile organic compounds (BVOCs) that are important precursors for tropospheric O 3 formation, affect methane lifetime, and also act as precursors for secondary organic aerosol (SOA).These latter interactions with SOA have a cooling effect (Arneth et al., 2010;Makkonen et al., 2012b;Paasonen et al., 2013).Condensation of monoterpenes (MTs), a group of BVOCs with large source strength from coniferous vegetation, on pre-existing particles increases the observed particle mass, as well as the number of particles large enough to act as cloud condensation nuclei (CCN; equivalent to particles > ca. 100 nm) at boreal forest sites (Tunved et al., 2006).For present-day conditions, Spracklen et al. (2008a) estimated a radiative cooling of −1.8 to −6.7 W per m −2 of boreal forest area from the BVOC-SOA interplay.
How future changes in MT emissions affect SOA growth and climate is very uncertain.This is partially because of the lack of process understanding of the various steps of aerosol formation and growth and interactions with cloud formation (Hallquist et al., 2009;Carslaw et al., 2010), and partially because the issue of how spatial patterns of changing emissions of atmospherically rapidly reactive substances translate into a changing patterns of radiative forcing, and then into a surface temperature change, has not yet been resolved (Shindell et al., 2008;Fiore et al., 2012).
The Russian boreal forest represents the largest continuous conifer region in the world.About one-third of this forested area (ca.730 × 10 6 ha) is dominated by larch (Shvidenko et al., 2007), in particular by the Larix gmelinii and L. cajanderii forests growing east of the Yenisei River on permafrost soils.Despite its vast expanse, the first seasonal measurements of MT emissions from Eastern Siberian larch have only recently been published (Kajos et al., 2013).Leaf MT emission capacities are highly species dependent; thus any model estimate of MT emissions from boreal larch forests that relies solely on generic BVOC emission parameterizations obtained from other conifer species will give inaccurate emission and related SOA aerosol number concentrations for this region (Spracklen et al., 2008a, b).We therefore provide here a first assessment of MT emission rates from the Eastern Siberian larch biome, combining measured emission capacities with a process-based dynamic vegetation model and quantitatively linking MT emissions and SOA formation.We use the observations and process models to assess climate change effects on future vegetation composition, BVOC emissions, and the concentration of particles of CCN size.We discuss how the climate impact of future SOA levels from changes in BVOC emissions across Eastern Siberia compares with changes in the regional CO 2 balance.The chief goal of the study is not to provide a full surface climate feedback quantification (for which today's global coupled modelling tools are insufficient) but rather to highlight the number of potentially opposing processes that need to be covered when doing so.

Site description, BVOC and aerosol measurements
Leaf BVOC emissions fluxes, above-canopy monoterpene concentration, and aerosol particle size and number concentration were measured during the growing season 2009 at the research station Spasskaya Pad, located ca. 40 km to the northeast of Yakutsk (62 • 15 18.4 N, 129 • 37 07.9 E) and centred in the Eastern Siberian larch biome (Kobak et al., 1996;Tchebakova et al., 2006).In the northern direction, no major pollution sources exist within hundreds of kilometres, the nearest mining areas are concentrated to the south and west of Yakutsk.The predominant air flow to the site is either from southern (via Yakutsk) or northern locations.Forest fires contribute to aerosol load in summer.
An eddy covariance tower for measurements of forestatmosphere exchange of CO 2 , water vapour, and sensible heat was established at Spasskaya Pad in the late 1990s (Ohta et al., 2001;Dolman et al., 2004) in a L. cajanderii forest growing on permafrost soil with an understory vegetation consisting of ericaceous shrubs.The forest has an average age of ca.185 years and canopy height is little less than 20 m.Maximum one-sided larch leaf area index (LAI) in summer is around 2 (Ohta et al., 2001;Takeshi et al., 2008).In 2009, leaf samples for BVOC analyses were taken, accessing the upper part of the canopy from a scaffolding tower located within a few hundred metres of the eddy flux tower (Kajos et al., 2013).Using a custom-made Teflon branch chamber, air filtered of O 3 was sampled onto Tenax-TA/Carbopack-B cartridges with a flow rate of 220 mL min −1 .A total of 5-12 samples were taken during the day from two trees on southfacing branches approximately 2 m below the tree top.The cartridge samples were stored at 5 • C during the campaigns, transported afterwards to Helsinki and thermally desorbed and analysed using a thermal desorption instrument (Perkin-Elmer TurboMatrix 650, Waltham, USA) attached to a gas chromatograph (Perkin-Elmer Clarus 600, Waltham, USA).For details on chamber, adsorbents, and laboratory measurements, see Haapanala et al. (2009), Ruuskanen et al. (2007), and Hakola et al. (2006).
Monoterpene concentrations and forest-atmosphere exchange fluxes were measured with a high-sensitivity quadrupole proton-transfer-reaction mass spectrometry (PTR-MS; Ionicon, Innsbruck, Austria) located in a hut at the foot of the eddy covariance tower.Sample air was drawn through a heated PFA tube using a 20 L min −1 flow from the inlet located at 30.3 m above ground.While reporting here on monoterpenes only, a range of masses, corresponding to BVOCs (e.g.isoprene, methanol, acetaldehyde), were sampled sequentially, with typical dwell times of 0.5 s and scanning sequences of around 4 s.Measurement set-up, disjunct eddy covariance flux calculations, and quality control followed Holst et al. (2010) could not be calibrated on-site.However, the instrument had been calibrated before and after the field campaign using a gas standard mixture from Ionimed (Innsbruck, Austria) using the same detector and instrument settings as during the field campaigns.Aerosol particles were continuously monitored with a scanning mobility particle sizer (SMPS) located at the foot of the eddy covariance tower, connected to a differential mobility analyser (DMA, medium Hauke type, custom built; for size segregation of aerosol particles) in front of a condensation particle counter (CPC, model 3010, TSI Inc., USA; for determining the number of the size segregated particles).The system was identical to the one described and evaluated in Svenningsson et al. (2008).Scans across the size range of 6-600 nm were completed every 5 min.The SMPS data were used to determine occasions of aerosol particle nucleation.The growth rates (GRs) were calculated from log-normal modes fitted to the measured particle size distribution following Hussein et al. (2005).The time evolution of the diameters at which the fitted modes peaked was inspected visually, and the GR was determined with linear least squares fitting to these peak diameters whenever a continuous increase in diameter was observed.In this analysis we calculated GRs for particles from 25 to 160 nm.
The source rate for condensing vapour (Q, molecules cm 3 s −1 ) was determined by calculating the concentration of condensable vapour needed to produce the observed GR (C GR , cm −3 ; Nieminen et al., 2010Nieminen et al., , 2014) ) and the condensation sink from the particle size distribution (CS, s −1 ; Kulmala et al., 2001).In steady state the sources and sinks for the condensing vapour are equal, and thus we determined the source rate as Q = C GR • CS.

Modelling of dynamic vegetation processes, permafrost, and BVOC emissions
We applied the dynamic global vegetation model (DGVM) LPJ-GUESS (Smith et al., 2001;Sitch et al., 2003), including algorithms to compute canopy BVOC emission following Niinemets et al. (1999), Arneth et al. (2007b), andSchurgers et al. (2009a) and permafrost as adopted from Wania et al. ( 2009) (Table 1).LPJ-GUESS simulates global and regional dynamics and composition of vegetation in response to changes in climate and atmospheric CO 2 concentration.Physiological processes like photosynthesis, autotrophic, and heterotrophic respiration are calculated explicitly; a set of carbon allocation rules determines plant growth.Plant establishment, growth, mortality, and decomposition, as well as their response to resource availability (light, water), modulate seasonal and successional population dynamics arising from a carbon allocation trade-off (Smith et al., 2001).Fire disturbance is included in the model (Thonicke et al., 2001).Similar to other DGVMs, a number of plant functional types (PFTs) are specified to represent the larger global vegetation units (Sitch et al., 2003).BVOC emissions models, whether these are linked to DGVMs or to a prescribed vegetation map, all rely on using emission potentials (E * , leaf emissions at standardized environmental conditions) or some derivatives in their algorithms.In LPJ-GUESS, production and emissions of leaf and canopy isoprene and monoterpenes are linked to their photosynthetic production, specifically the electron transport rate, and the requirements for energy and redox equivalents to produce a unit of isoprene from triose phosphates (Niinemets et al., 1999;Arneth et al., 2007b;Schurgers et al., 2009a).A specified fraction of absorbed electrons used for isoprene (monoterpene) production (ε) provides the link to PFT-specific E * (Arneth et al., 2007a); in case of monoterpenes emitted from storage, an additional correction is applied to account for their light-dependent production (taking place over parts of the day) and temperature-driven emissions (taking place the entire day) (Schurgers et al., 2009a).Half of the produced monoterpenes were stored, whereas the other half was emitted directly (Schurgers et al., 2009a).Values for E * were similar to the global parameterization for most of the model's PFTs (Schurgers et al., 2009a), with the exception of boreal needle-leaf summergreen (BNS) "larch" PFT (see results).
Contrasting the stimulation of BVOC emissions in a warmer and more productive environment, higher CO 2 concentrations have been shown to inhibit leaf isoprene production.Even though the underlying metabolic mechanism is not yet fully understood, this effect has been observed in a number of studies (for an overview see Fig. 6 in Arneth et al., 2011).Due to limited experimental evidence, whether or not a similar response occurs in monoterpene producing species cannot yet be confirmed, especially in species that emit from storage.The model is set up to test this hypothesis.Arneth et al. (2007a) proposed an empirical function for CO 2 inhibition, based on the ratio of leaf internal CO 2 concentration at a standard atmospheric CO 2 level (taken as 370 ppm) and at the given atmospheric CO 2 levels of the simulation year (both calculated for non-water-stressed conditions); the relationships have been shown since to fit an updated compilation of observations well (Arneth et al., 2011).The algorithm that describes the CO 2 inhibition of BVOC emission can either be enabled or disabled in the model (Arneth et al., 2007a) and simulations results thus compared (see Fig. A1, Appendix; Table 2).
LPJ-GUESS was recently expanded with a permafrost module following Wania et al. (2009;Miller and Smith, 2012) in which a numerical solution of the heat diffusion equation was introduced.The soil column in LPJ-GUESS now consists of a snow layer of variable thickness, a litter layer of fixed thickness (5 cm), and a soil column of depth 2 m (with sublayers of thickness 0.1 m) from which plants can extract non-frozen water above the wilting point.A "padding" column of depth 48 m (with thicker sublayers)  N) and represent an area of 1.2 × 10 7 km 2 .NPP global (given as a reference value) is global vegetation net primary productivity (Pg C a −1 ).BVOCs are in Tg C a −1 , CO 2 -C fluxes in Pg C a −1 ; C pools are in PgC.Simulations for monoterpene emissions for the boreal needleleaf summergreen (BNS) plant functional type compared three cases (indicated as different subscripts for "Total_MT BNS "), using maximum (9.6 µgC g −1 h −1 ) and minimum (1.9 µgC g −1 h −1 ) values for E * measured in Spasskaya Pad (see text) and E * = 6.2 µgC g −1 h −1 as a weighted average from all observations at the Spasskaya Pad location.For BVOCs, CO 2 inhibition was switched on and off (Arneth et al., 2007b(Arneth et al., ). 1981(Arneth et al., -2000(Arneth et al., 2031(Arneth et al., -2050  is also present beneath these three layers to aid in the accurate simulation of temperatures in the overlying compartments (Wania et al., 2009).Soil temperatures throughout the soil column are calculated daily, in addition to change in re-sponse to changing surface air temperature and precipitation input and the insulating effects of the snow layer and phase changes in the soil's water.
Here we run the model with 0.5 • spatial resolution, using climate and atmospheric CO 2 as driving variables as described in the literature (Smith et al., 2001).Values describing growth and survival of the BNS (larch) PFT were adopted from previous studies (Sitch et al., 2003;Hickler et al., 2012;Miller and Smith, 2012), but with the degree-day cumulative temperature requirements on a 5 • basis (GDD5) to attain full leaf cover reduced from 200 to 100 (Moser et al., 2012).Minimum GDD5 to allow establishment was set to 350 resulting in establishment of seedlings in very cold locations.Soil thermal conductivity was 2 W m −1 K −1 .The model was spun up for 500 years to 1900 values using CO 2 concentration from the year 1900 and repeating de-trended climate from 1901 to 1930 from CRU (Mitchell and Jones, 2005).Historical (20th century) simulations used observed CO 2 concentrations and based on variable CRU climate.Simulations for the 21st century were based on ECHAM climate, using RCP 8.5 emissions (Riahi et al., 2007).The model requires daily radiation, precipitation, and maximum and minimum air temperatures as input (Arneth et al., 2007b).Climate output from ECHAM was interpolated to the CRU half-degree grid, and monthly values interpolated to daily ones (see Ahlström et al., 2012, and references therein).These daily fields were then bias corrected using the years 1961-1990 as reference period, as in Ahlström et al. (2012).CO 2 inhibition of BVOC emissions were switched on and off in separate simulations to assess the sensitivity of our results to this process.Totals across Siberia were calculated for a grid box that ranged from 46 to 71 • N and 76 to 164 • E. Simulated changes in total carbon uptake or losses were translated into radiative forcing following (IPCC, 2007), assuming a 50 % uptake in oceans in case of a net loss to the atmosphere (Sitch et al., 2007).

Modelling aerosols and CCN
To model the effect of BVOCs on CCN concentrations, we use the global aerosol-climate model ECHAM5.5-HAM2(Zhang et al., 2012).ECHAM5.5-HAM2includes the aerosol components of black carbon, organic carbon, dust, sea salt, and sulfate (Table 1) and describes the aerosol size distribution with seven log-normal modes.The microphysics module M7 (Vignati et al., 2004) includes nucleation, coagulation, and condensation.In this study, we use the ECHAM5.5-HAM2version with activation type as described in Makkonen et al. (2012a, b).For simulating secondary organic aerosol, we use the recently developed SOA module (Jokinen et al., 2015).The SOA module explicitly accounts for gas-phase formation of extremely low-volatility organic compounds (ELVOCs) from monoterpene oxidation.The module implements a hybrid mechanism for SOA formation: ELVOCs are assumed to condense to the aerosol population according to the Fuchs-corrected condensation sink, while semi-volatile organic compounds (SVOCs) are partitioned according to organic aerosol mass.While simulated ELVOCs are able to partition more effectively to nucle-ation and Aitken mode, hence providing growth for nucleated particles to CCN size, SVOCs primarily add organic mass to accumulation and coarse aerosol modes.A total SOA yield of 15 % from monoterpenes is assumed (Dentener et al., 2006).While similar assumption on total SOA yield is applied by most aerosol-climate models, the simulated SOA is likely to be underestimated (e.g.Tsigaridis et al., 2014).
ECHAM5.5-HAM2 was run with different BVOC emission scenarios that were simulated for the years 2000 and 2100 off line with LPJ-GUES (see previous section).The model uses T63 spectral resolution with 31 vertical hybrid sigma levels.The simulations apply present-day oxidant fields as in Stier et al. (2005).All simulations are initiated with a 6-month spin-up, followed by 5 years of simulation for analysis.The model climate is nudged towards ERA-40 reanalysis year 2000 meteorology, an approach that is widely used in aerosol-climate assessments (K.Zhang et al., 2014).Present-day wildfire and anthropogenic aerosol and precursor emissions are applied for all simulations (Dentener et al., 2006).One of the foci here is a BVOC, comparing presentday and future BVOC emissions (choosing a conservative estimate of E * = 1.9 µg C m −2 (leaf) h −1 ; see results section for further details on E*) but keeping other emissions constant (oxidant fields and nudging meteorology are same for both 2000 and 2100).The emissions of dust and sea salt are modelled interactively (Zhang et al., 2012).
The analysis of model results includes total particle number concentrations and CCN at 1 % supersaturation (CCN(1 %)), since it reflects the changes in Aitken mode concentrations and local changes in precursor emissions.While "realistic" supersaturations are generally lower, choosing CCN(1.0 %) concentration provides the upper limit for CCN concentrations.The simulations are also used to assess the radiative effects of SOA.In the simulations, the aerosol concentrations are interactively coupled to the cloudmicrophysics scheme (Lohmann et al., 2007) and to the direct aerosol radiative calculation.The aerosol indirect effect is evaluated as a change in cloud radiative forcing ( CRF).The direct aerosol effect accounts only for clear-sky shortwave forcing ( CSDRF).The radiative effects are calculated as differences from two time-averaged 5-year simulations as Subscripts "2000" and "2100" denote that BVOC emissions from this year were used, while other model conditions were based on present-day values.

Present-day expanse of larch forest and BVOC emissions
LPJ-GUESS reproduces the present-day circumpolar permafrost distribution (Fig. 1; shown as circumpolar map for comparison with Tarnocai et al., 2009) and, with the exception of the Kamchatka peninsula, simulates also the expanse of the larch-dominated forests in Eastern Siberia (Fig. 1; Miller and Smith, 2012;Wagner, 1997).Maximum LAI calculated by the model for the Spasskaya Pad forest (62 • 15 18.4 N, 129 • 37 07.9 E; 220 m a.s.l.), where the BVOC measurements were obtained, was 2.0 (averaged over years 1981-2000; not shown) and is in good agreement with the measured values during that period (1.6; Takeshi et al., 2008).Total present-day modelled soil C pools over the top 2 m in Eastern Siberia are 216 Gt C and, for circumpolar soils summed for latitudes above 40 • N, 454 Gt C (Table 2).Recent estimates of C stored in northern latitude soils affected by permafrost were 191, 495, and 1024 Gt C in the 0-30, 0-100, and 0-300 cm soil layer respectively, based on extrapolating observations stored in the Northern Circumpolar Soil Carbon Database (Tarnocai et al., 2009).These numbers indicate that the values calculated with LPJ-GUESS are lower than observation-based ones, most likely underestimating C density in particular in the soil layers below a few tenths of centimetres.Kajos et al. (2013) measured for the first time MT E * from L. cajanderii.Their measurements, taken over an entire growing season at Spasskaya Pad, suggested values of E * ranging from 1.9 µg C m −2 (leaf) h −1 at the lower end to 9.6 µg C m −2 (leaf) h −1 at the upper.Applying a weighted measured-average E * of 6.2 µg C m −2 (leaf) h −1 in LPJ-GUESS led to average summer daily monoterpene emissions of 2.9 (± 0.8, 1 standard deviation, June) mg C m −2 d −1 , and 2.2 (± 0.8, July) mg C m −2 d −1 for the grid location representing the Spasskaya Pad site, for the same year of measurements as reported in Kajos et al. (2013).These values are within 25 and 10 % of measured values (3.3 ± 2.9 mg C m −2 d −1 , June; 2.4 ± 1.6 mg C m −2 d −1 , July), even though the modelled day-to-day variation was smaller, which is expected when applying grid-averaged climate as model input.By comparison, for a boreal Scots pine forest in southern Finland the average June and July monoterpene emissions were somewhat larger than the values for larch (3.8 and 5.1 mg C m −2 d −1 respectively; Rantala et al., 2015).For a Larix kaempferi-dominated temperate forest in Japan, Mochizuki et al. (2014) extrapolated, based on their measurements, summertime maxima of ca. 10 mg C m −2 d −1 .Across the Siberian larch biome, applying E * of 6.2 µg C m −2 increased simulated total present-day MT emissions from 0.11 TgCa −1 (as in Schurgers et al., 2009a) to 0.21 TgCa −1 , or to 0.42 TgCa −1 when the maximum E * was used (Table 2).The observed range in E * , and the calculated range in total emissions across Siberia, might reflect variability in tree microclimate or genetic variability or might have been induced by (undetected) mechanic or biotic stress during the time of measurements (Staudt et al., 2001;Bäck et al., 2012;Kajos et al., 2013).While our data are insufficient to make a finite suggestion of L. cajanderi E * , the measurements provide evidence for potentially substantially higher MT emissions from Siberian larch than previous estimates.

Present-day aerosols and links to BVOCs
New particle formation events (Fig. 2a) were observed regularly at Spasskaya Pad.The calculated volumetric source rates of condensing vapours (Q), the product of vapour concentration required for the observed particle growth rate and particle loss rate (Kulmala et al., 2005), increased exponentially with temperature (Fig. 2b).MT concentrations increased with temperature as well, with a slope relatively similar to that found for the Q vs. T relationship (Fig. 2c).Consequently, a positive relationship emerged between Q and MT concentration (Fig. 2d), which supports previous field and laboratory evidence that MTs and their oxidation products are a main precursor to the observed particle formation and growth.
Figure 2d shows the connection between the BVOC concentration and the formation rate of vapours causing the growth of the aerosol particles.Even though the monoterpene concentrations were measured above and the aerosol growth rates below the canopy, the observed correlation indicates that BVOC concentration is an important contributor to the regional aerosol growth and supports the theory that the condensation of organic vapour is largely responsible for the formation of secondary organic aerosol (Hallquist et al., 2009;Carslaw et al., 2010).Substantial within-canopy chemical reactions would be expected to worsen the relationship.The correlation depicted in Fig. 2d is determined in particular by the formation of secondary organic aerosol on preexisting aerosol particles, whereas the nucleation rate of new aerosol particles seems not to be dominated by the landscapescale emissions and surface concentrations of BVOCs.For instance, most nucleation events in a Scots pine dominated landscape in Finland have been found in spring, when measured monoterpene concentrations in the near-surface were about one-tenth of the summertime maximum (ca.60 ppt, vs. up to 500 ppt; Haapanala et al., 2007;Lappalainen et al., 2009).In our study, we found MT concentrations of similar magnitude to these.
In contrast to temperature and BVOC concentrations, levels of radiation, which can be considered a surrogate for the concentration of the OH radical (OH q ), did not affect Q (Fig. 2b), even though OH q has been considered an important player in aerosol formation.Rohrer and Berresheim (2006) showed a strong correlation between solar ultraviolet radiation and OH q concentration at the Hohenpeissenberg site in Germany.Furthermore, Hens et al. (2014) demonstrated that the daytime OH q concentrations in (especially) boreal forest depend on solar radiation, both above and below the canopy.Hence, the poor relation between the source rate of condensing vapour and orders-of-magnitude variation in levels of radiation (Fig. 2b) indicates that OH-radical concentration did not have a major impact on Q.This agrees with the findings by Ehn et al. (2014) that ozone instead of OH q is an important, if not the main, atmospheric agent oxidizing organic vapours into a chemical form that condenses on particle surfaces.Since the relative variation in ozone concentrations is much smaller than in BVOC (or OH) concentrations (Hens et al., 2014), the similarity in the dependencies of Q and MT concentration and temperature (Fig. 2b and c) is in favour of a more significant role of ozone than of OH in the formation of condensable vapours.In general, our results indicate that factors and processes besides the concentrations of SO 2 and OH q seem to limit aerosol production in unpolluted environments (Kulmala et al., 2005).

Future carbon pools, vegetation distribution, and BVOC emissions in Siberia
In a warmer environment with higher atmospheric CO 2 levels, the simulations indicated a drastically reduced area of permafrost in Siberia (Fig. 1).Total net primary productivity in the simulated domain increased from an annual average of 3.5 to 5.9 PgC a −1 at the end of the 21st century.An overall C loss of 100 PgC assumed to be in the form of CO 2 (since the model does not yet include a dynamic surface hydrology which would be necessary to assess changing methane emissions) at the end of the 21st century, compared to presentday conditions, was calculated from the shrinking Siberian areas of permafrost (Table 2).However, warming and higher levels of atmospheric CO 2 led also to increasing LAI and to larch-dominated areas showing the expected north and north-eastwards shift (Fig. 1) compared to present-day climate (Miller and Smith, 2012).The carbon uptake in expanding vegetation into permafrost-free areas, combined with enhanced productivity across the simulation domain overcompensates for the losses from C pools in permafrost areas (Table 2).Without CO 2 inhibition of BVOC emissions future MT emissions were, as expected, notably enhanced: directly as a result of warmer leaves, but augmented by the future higher LAI of larch and evergreen conifers (Fig. 1d; Table 2, Fig. A1).Since the emissions scale with the emission factors applied, the proportional increase between present-day and future climate conditions is independent of the value of E * .Whether or not leaf MT emissions are inhibited by increasing atmospheric CO 2 levels to a similar degree to what was found for isoprene is difficult to assess from today's lim- ited number of studies (e.g.Niinemets et al., 2010, and references therein).We included both simulation results in Table 2 since similarities in the leaf metabolic pathways of isoprene and MT production suggest such an inhibition, but possibly this effect does not become apparent in plant species where produced MTs are stored unless the storage pools become measurably depleted by the reduced production.By contrast, species emitting MTs in an "isoprene-like" fashion immediately after production should more directly reflect CO 2 inhibition.Evergreen conifers typically emit most MTs from storage pools, although recent experiments have shown that some light-dependent emissions also contribute to total emission fluxes.Accordingly, based on the leaf-level measurements, larch could follow a hybrid pattern between emission after production and from storage (Kajos et al., 2013).Without accounting for CO 2 inhibition, MT emissions across the model domain more than doubled (Fig. 1; Table 2) by 2100, as a consequence of higher emissions per leaf area due to warmer temperatures and of the larger emitting leaf area in response to higher photosynthesis.With CO 2 inhibition included, simulated changes were negligible, similar to what was shown in previous BVOC simulations with the model (Arneth et al., 2007a(Arneth et al., , 2008)).

Discussion
Boreal vegetation has been shown to respond to the recent decades' warming and increasing atmospheric CO 2 levels with a prolonged growing season and higher maximum LAI, similar to patterns in our simulations (Piao et al., 2006).The calculated enhanced biomass growth is in line with experimental evidence of higher C in plant biomass in warming plots at tundra field sites (Elmendorf et al., 2012;Sistla et al., 2013).In Siberian mountain regions, an upward movement of vegetation zones has been recorded already (Soja et al., 2007), while the analysis of evergreen coniferous undergrowth abundance and age shows spread of evergreen species, especially Pinus siberia, into Siberian larch forest (Kharuk et al., 2007).These observations thus support the modelled shift in vegetation zones and the change in vegetation type composition and productivity.Likewise, other models with dynamic vegetation also have shown a strong expansion of broadleaved forests at the southern edge of the Siberian region in response to warming (Shuman et al., 2015).
Warming and thawing of permafrost soils is being observed at global monitoring network sites, including in Russia (Romanovsky et al., 2010).Estimates of carbon losses from northern wetland and permafrost soils in response to 21st century warming range from a few tens to a few hundreds Pg C, depending on whether processes linked to microbial heat production, thermokarst formation, and surface hydrology, winter snow cover insulation, dynamic vegetation, C-N interactions, or fire are considered (Khvorostyanov et al., 2008;Schuur et al., 2009;Arneth et al., 2010;Koven et al., 2011;Schneider von Deimling et al., 2012).For instance, a modelled range of 0.07-0.23W m −2 forcing associated with a 33-114 Pg CO 2 -C loss from permafrost regions was found for a simulation study that was based on the RCP 8.5 climate and CO 2 scenarios, but excluding full treatment of vegetation dynamics (Schneider von Deimling et al., 2012).In a recent literature review, Schaefer et al. (2014) found a range from cumulative 46 to 435 CO 2 equivalents (accounting for CO 2 and CH 4 ), or 120 ± 85 GtC by 2100, in response to different future warming scenarios and modelling approaches.In our simulation, the CO 2 -C loss from the decreasing Siberian permafrost region would be equivalent to a 0.13 additional W m −2 forcing in 2100 (see methods).Likely, this number is too low since the model does not include thermokarst processes, which can facilitate rapid thaw (Schaefer et al., 2014, and references therein).The modelled carbon loss was offset when taking into account vegetation dynamics and processes across the entire Siberian study domain (Table 2), including a shift in PFT composition, and enhanced productivity especially in the southern regions, such that the overall carbon uptake including enhanced net primary productivity and expanding woody vegetation resulted in a small negative (−0.09 W m −2 ) effect.
LPJ-GUESS is a second-generation DGVM (Fisher et al., 2010) and includes plant demography, such that forest successional dynamics and competition for water and light between individual age cohorts are treated explicitly (Smith et al., 2001).The forest growth dynamics thus differentiate between early successional, short-lived species that invest in rapid growth and shade-tolerant trees with resource allocation aimed towards longer-lived growth strategies.As a result, the model's PFTs can be mapped to tree species when required information for model parameterization is available.
The process-based treatment of resource competition such as for light and water has been shown to lead to a realistic growth response and distribution under present-day climate conditions (Arneth et al., 2008;Schurgers er al., 2009b), which should also hold in future and past climates (Miller et al., 2008;Schurgers et al., 2009b).This feature also provides a distinct advantage when applying the necessary BVOC emission capacities that are based on species (rather than functional-type) average values (Arneth et al., 2008;Schurgers et al., 2009b;Niinemets et al., 2010).In earlier simulations (Schurgers et al., 2009a), a generic emission potential of E * = 2.4 µg C m −2 (leaf) h −1 was adopted for the BNS PFT based on a recommendation in Guenther et al. (1995), which at that time did not include observations from any larch species.Here we demonstrate not only that a range of measured larch E * (see Table 2) introduces large uncertainty in total MT emissions from Siberia but also that it is fundamental to apply dynamic vegetation growth response (rather than static maps) for BVOC emissions estimates in changing environments.
Monoterpene compounds can be emitted either directly following their synthesis in the chloroplast, in an "isoprenelike" fashion, or from storage pools, resulting in an emission pattern that is independent of light availability.The observed emissions of monoterpenes by larch possibly reflect a hybrid pattern between emission directly after synthesis in the chloroplast and emission from storage pools, as has also been found for other coniferous species (Schurgers et al., 2009a).The needle-level measurements by Kajos et al. (2013) on larch indicated a combined light and temperature response, even though a robust differentiation to a temperature-only model was not possible due to the limited sample size.An earlier study by Ruuskanen et al. (2007) on a 5-year-old L. sibirica tree indicated a better performance of the temperature-only emission model for monoterpene species compared to the light and temperature approach.
Multiple interacting processes could alter monoterpene emissions in the future.Irrespective of the relative roles of light vs. temperature dependence, a change in MT concentrations and hence partial pressure of MT in stored pools, for instance in response to long-term warming, would affect emission capacities.Changes in measured E* when investigated over the course of a growing season have been reported and could be related to a changing production rate (Niinemets et al., 2010).Likewise, observed profiles of E * within tree canopies appear related to not only changes in leaf area-to-weight ratios along the canopy light and temperature gradients but also to varying production rates (Niinemets et al., 2010).Emission capacities in Q. ilex leaves adapted to warm growth environment were notably enhanced (Staudt et al., 2003), but the experimental basis for an acclimation response of BVOC emissions to temperature remains remarkably poor (Penuelas and Staudt, 2010) and is indicative of the general lack of global modelling studies accounting for possibly acclimation of process responses to environmental changes (Arneth et al., 2012).In our simulations we aim to provide a range of a possible plastic BVOC-CO 2 responses by switching the direct CO 2 inhibition on and off for both isoprene and monoterpene, but we do not account for other acclimation processes.
The assessment of climate effects of changes in the CO 2 -C balance vs. those of BVOC-SOA interactions is challenging, since the translation of regional changes in emissions of atmospherically reactive species into related radiative forcing and then into a response in the climate system is highly nonlinear and poorly understood (Shindell et al., 2008;Fiore et al., 2012) sonen et al., 2013), but a simple extrapolation based on the region's growing season temperature increases of ca.5.5 K simulated at the end of the 21st century in our study with the ECHAM GCM does not account for the important nonlinearities in the system.Present-day CCN (1.0 %) concentrations over Siberia were estimated to vary from extremely low values of less than 50 cm −3 north of 60 • N to a few hundred per cc in the southern part of Siberian domain (Fig. 3).In order to put measurements and model simulations into context, simulated CCN concentrations (at the Spasskaya Pad location) were evaluated against observations during May-August, using particle diameter (D p ) > 100 nm as proxy for CCN, since CCN at different supersaturations were unavailable in the observations.The model reproduces the observed May-August average CCN concentration and CCN maximum location in July (not shown), but the seasonal variation was overestimated in the simulations.ECHAM-HAM indicates a transition from very clean spring aerosol population of ∼ 100 cm −3 to high July concentrations ranging from 800 to 1200 cm −3 in the Yakutsk region.By contrast, obser-vations show only moderate monthly CCN variability from 550 cm −3 in May to 750 cm −3 in July.While the simulated low spring concentrations likely reflect unaccounted-for anthropogenic emissions, the simulated high summer concentrations result from strong wildfire emissions in the region in the applied emission inventory (see below).
Whether or not BVOCs can increase the availability of CCN depends on the availability of sub-CCN-sized particles (O'Donnell et al., 2011).In the future, a scenario of decreasing anthropogenic emissions could lead to a strong decrease in calculated atmospheric SO 2 concentrations and also in particle nucleation (Makkonen et al., 2012a).In the model experiments anthropogenic primary emissions are introduced as 60 nm particles; hence condensation of sulfuric acid and organic vapours is generally needed in order to grow these particles to CCN sizes.However, the modelled primary particle emissions are dominated by wildfires, which are assumed to inject large particles with 150 nm diameter.SOA formation only partly enhances the survival of small particles by providing additional growth (Makkonen et al., 2012a), but partly also suppresses it by increasing the coagulation sink for small particles (Fig. A2, lower left panel; see also O'Donnell et al., 2011).
The assumption of unchanging oxidant fields induces some uncertainty for future simulations and inconsistency with present-day simulations with varying biogenic emissions, since both anthropogenic and biogenic emissions are likely to modify the atmospheric oxidative capacity.Nudging towards reanalysis meteorology establishes evaluation of BVOC-aerosol coupling with unchanged meteorological fields but restricts the model in terms of aerosol-climate feedbacks, since e.g.nudging future climate simulations with present-day meteorological winds is based on the assumption that wind direction and speed, etc., are not changing.
When only BVOC emissions were changed between the present day and levels simulated for climate in 2100, the relatively higher emission of BVOCs leads to substantially increased aerosol growth rates over a large part of the The areas are averaged over Siberia, and the BVOC emissions for years 2000 to 2100 (example is forE * = 1.9).Areas were separated by wildfire emissions (using an emission limit of 10 −11 kg m −2 s −1 ).In the Siberian domain, accumulation mode includes over 85 % of organic aerosol, and the absolute changes in SOA are also dominated by accumulation mode.However, the SOA condensation increase until year 2100 is essential for nucleation and Aitken mode growth.
Siberian domain.This was the case even though we chose the conservative estimate based on the low measured E * of 1.9 µg C m −2 h −1 .However, GR is not the only factor determining levels of CCN.Increased aerosol mass due to increased SOA formation led also to an increase in the condensation sink and eventually to even decreased particle formation rates in some regions (Fig. A2, lower right panel).These competing effects of increased growth and increased sink are essential for quantifying the importance of the cloud albedo forcing feedback.We can also show that the patterns of changes in CCN in response to future BVOC emissions are additionally affected by changes in the aerosol background, which strongly influences the indirect aerosol effect of SOA.
In large parts of Siberia, the simulated BVOC oxidation products condense on CCN-sized aerosols already present from wildfires.When simulation results were separated into regions of low and high wildfire emissions (Fig. 4), areas of low wildfire activity had a relatively large increase in SOA formation (60 %) in nucleation mode (D p < 10 nm).The relative increases in SOA formation in Aitken, accumulation, and coarse modes were 50, 31, and 40 % respectively.However, the distribution of BVOC oxidation products was rather different in areas of high wildfire activity: the condensation of SOA depends on surface area and organic mass of the population, both of which are shifted towards larger modes in wildfire-intensive areas.SOA formation in coarse mode was more than doubled, while SOA in nucleation mode decreased by 30 % due to a decrease in nucleation rates and increase in vapour sink in large aerosol modes.
It is clear that the effect of increased BVOC emission on particle population has distinct effects depending on existing background aerosol distribution.Moreover, CCN at 1.0 % supersaturation was used even though "realistic" supersaturations are generally lower.The CCN(1.0 %) concentration therefore provides an upper limit for CCN concentration.
In the aerosol model, neither the simulated CCN(1.0 %) nor CCN(0.2 %) corresponds clearly to either Aitken or accumulation modes.CCN at 0.2 % would reflect larger aerosols, and hence the changes in CCN(0.2 %) would be less sensitive to aerosol and precursor sources (see corresponding Fig. A3).
Averaged over Siberian areas of low wildfire activity, the median (mean) increase of CCN(0.2 %) was calculated to be 1 % (7 %) due to BVOC emissions changes from present-day levels to the end of the 21st century, while areas of high wildfire emission lead to median (mean) increase of 0.3 % (0.5 %).
Even though the Siberian MT emissions more than double until 2100 (Table 2), the increasing wildfire emissions and decreasing new particle formation due to reductions in anthropogenic SO 2 largely offset the effect of increased BVOC emissions on CCN concentration.In wildfire plumes, the simulated CCN concentrations were high even without BVOC-induced growth of smaller particles.The radiative effect due to BVOC emission change between ca.2000 and ca.2100 was estimated from ECHAM-HAM simulations averaged over 5 years.The increase in BVOC emissions leading to additional secondary organic aerosol induces a −0.2 W m −2 change in direct clear-sky aerosol forcing over the Siberian domain at the end of the 21st century.Furthermore, the increase in CCN concentrations leads to a strengthening of the cloud radiative effect by −0.5 W m −2 (Table 3).These changes in radiative fluxes only take into account the changing BVOC emission, and the potential concurrent changes in anthropogenic and wildfire emissions might decrease the simulated radiative effect of biogenic SOA (Carslaw et al., 2013).

Implications, limitations, and future progress
Up to now, studies that investigate the role of terrestrial vegetation dynamics and carbon cycle in the climate system typically account solely for CO 2 , while studies that look at BVOC-climate interactions often ignore other processes, especially interactions with vegetation dynamics or the CO 2 balance of ecosystems.However, for understanding the full range of interactions between atmospheric composition, climate change, and terrestrial processes we need a much more integrative perspective.Our analysis seeks to provide an example of how to quantify a number of climatically relevant ecosystem processes in the large Eastern Siberian region in a consistent observational and modelling framework that accounts for the multiple interactions between emissions, vegetation, and soils.It poses a challenge to combine effects of well-mixed greenhouse gases and locally constrained, shortlived substances.On a global scale, the opposing estimates in radiative effects from ecosystem-CO 2 and BVOC-SOA interactions are minuscule but it is to be expected that some of the forcing effects from SOA could lead to a notable change in regional temperatures.Clearly our numbers are uncertain, but they pinpoint the necessity for assessing surfaceatmosphere exchange processes comprehensively in climate feedback analyses.
While doing so, we are aware of the fact that a number of additional processes are not included in our analysis.For instance, it remains to be investigated whether a similar picture would emerge when additional feedback mechanisms are taken into consideration, e.g.SOA formation from isoprene (Henze and Seinfeld, 2006) or effects of atmospheric water vapour on reaction rates and aerosol loads, or that some of the SOA might like to partition more to the gas phase in a warmer climate.Likewise, neither the albedo effect of northwards migrating vegetation (Betts, 2000;W. Zhang et al., 2014), nor changes in the hydrology (which affect CH 4 and N 2 O vs. CO 2 fluxes), nor changes in C-N interactions (Zaehle et al., 2010) are considered here, which would require a coupled Earth system model that combines a broad range of dynamically varying ecosystem processes with full treatment of air chemistry and aerosol interactions.Quantifying the full range of terrestrial climate feedbacks, either globally or regionally, with consistent model frameworks that account for the manifold interactions is not yet possible with today's modelling tools.

Figure 1 .
Figure 1.Simulated maximum summer monthly leaf area index (LAI; a, b) and July emissions of monoterpenes (c, d; mgC m −2 month −1 ) from Eastern Siberian larch.The latter were calculated applying emission factors of 6.2, obtained from the measurements at Spasskaya Pad.(e, f) Maximum permafrost thaw depth (August), shown here as the circumpolar map for comparison with Tarnocai et al. (2009).Values are averages for a simulation 1981-2000 (a, c, e) and for 2081-2100 (b, d, f), applying climate and CO 2 concentrations from ECHAM-RCP8.5.Emissions in (c, d) do not account for direct CO 2 inhibition (see also Fig. A1).

Figure 2 .
Figure 2. Particle growth rates obtained from particle number size distribution (a, example from 10 June 2009).The colours indicate the measured concentrations (dN/d log D p , cm −3 ) of particles with different diameters (D p , nm) over the course of a day, small circles are mean diameters of concentration modes fitted for each measurement, and the temporal change of these diameters is represented with black lines from which the growth rate is calculated.Panel (b) shows the calculated volumetric source rates of condensing vapours (Q, molecules cm −3 s −1 ; 10 min resolution) as a function of air temperature ( • C) for all identified growth periods (one data point is obtained for each fitted growth rate, e.g. from (a) three data points would have been extracted); data are separated by levels of photosynthetically active radiation (PAR).(c) Monoterpene concentrations (half-hourly data) measured above the canopy vs. temperature measured at the same level (data separated by PAR, the data points overlapping with determined growth rate in b are indicated by encircled symbols), and relationship between volumetric source rate of condensing vapours and monoterpene concentration (d; data separated by particle diameter).Data points in (d) correspond directly to encircled symbols in (b).

Figure 3 .
Figure 3. Annual average boundary-layer CCN (1.0 %) concentration (cm −3 ) in Siberia with present-day anthropogenic and BVOCs (for BNS: E * = 1.9) emissions (left panel), and changes in CCN (1.0 %; right panel) concentration due to increase in BVOC emission between years 2000 and 2100 (simulations with CO 2 inhibition off).Areas with statistical significant changes in CCN are indicated by dots.The statistical analysis is based on monthly average CCN concentrations from 5 years of simulated data, and statistical significance of the CCN anomaly is evaluated using a two-sample t test, without assuming equal variance between the two populations.
. Based on a synthesis of measured aerosol number number concentrations and size distribution combined with boundary-layer growth modelling, Paasonen et al. (2013) estimated a growing-season indirect radiative cloud albedo feedback of −0.5 W m −2 K −1 for the Siberian larch region.The observation-based indirect feedback factors exceeded direct ones by roughly an order of magnitude (Paa

Figure 4 .
Figure 4. Relative increase in SOA mass, simulated by ECHAM5-HAM in different aerosol size modes due to BVOC emissions increase from the year 2000 to 2100.The areas are averaged over Siberia, and the BVOC emissions for years 2000 to 2100 (example is forE * = 1.9).Areas were separated by wildfire emissions (using an emission limit of 10 −11 kg m −2 s −1 ).In the Siberian domain, accumulation mode includes over 85 % of organic aerosol, and the absolute changes in SOA are also dominated by accumulation mode.However, the SOA condensation increase until year 2100 is essential for nucleation and Aitken mode growth.

Table 1 .
Overview of modelled processes and model-specific features.For further details see text.

Table 2 .
Simulated changes in net primary productivity, BVOC emissions, and C pool size in vegetation and soils.Unless stated otherwise, values are for the simulated Siberian domain (76-164 • E, 46-71

Table 3 .
Simulated changes in radiative effects due to change in BVOC emission between years 2000 and 2100, averaged over Siberian domain, Northern Hemisphere, and globally.CRF is cloud radiative forcing; CSDRF is direct aerosol effect, which accounts only for clear-sky short-wave forcing.