Challenges in modelling isoprene and monoterpene emission dynamics of Arctic plants: a case study from a subarctic tundra heath

Abstract. The Arctic is warming at twice the global average speed, and the warming-induced increases in biogenic volatile organic compounds (BVOCs) emissions from Arctic plants are expected to be drastic. The current global models' estimations of minimal BVOC emissions from the Arctic are based on very few observations and have been challenged increasingly by field data. This study applied a dynamic ecosystem model, LPJ-GUESS, as a platform to investigate short-term and long-term BVOC emission responses to Arctic climate warming. Field observations in a subarctic tundra heath with long-term (13-year) warming treatments were extensively used for parameterizing and evaluating BVOC-related processes (photosynthesis, emission responses to temperature and vegetation composition). We propose an adjusted temperature (T) response curve for Arctic plants with much stronger T sensitivity than the commonly used algorithms for large-scale modelling. The simulated emission responses to 2 °C warming between the adjusted and original T response curves were evaluated against the observed warming responses (WRs) at short-term scales. Moreover, the model responses to warming by 4 and 8 °C were also investigated as a sensitivity test. The model showed reasonable agreement to the observed vegetation CO2 fluxes in the main growing season as well as day-to-day variability of isoprene and monoterpene emissions. The observed relatively high WRs were better captured by the adjusted T response curve than by the common one. During 1999–2012, the modelled annual mean isoprene and monoterpene emissions were 20 and 8 mg C m−2 yr−1, with an increase by 55 and 57 % for 2 °C summertime warming, respectively. Warming by 4 and 8 °C for the same period further elevated isoprene emission for all years, but the impacts on monoterpene emissions levelled off during the last few years. At hour-day scale, the WRs seem to be strongly impacted by canopy air T, while at the day–year scale, the WRs are a combined effect of plant functional type (PFT) dynamics and instantaneous BVOC responses to warming. The identified challenges in estimating Arctic BVOC emissions are (1) correct leaf T estimation, (2) PFT parameterization accounting for plant emission features as well as physiological responses to warming, and (3) representation of long-term vegetation changes in the past and the future.

Abstract. The Arctic is warming at twice the global average speed, and the warming-induced increases in biogenic volatile organic compounds (BVOCs) emissions from Arctic plants are expected to be drastic. The current global models' estimations of minimal BVOC emissions from the Arctic are based on very few observations and have been challenged increasingly by field data. This study applied a dynamic ecosystem model, LPJ-GUESS, as a platform to investigate short-term and long-term BVOC emission responses to Arctic climate warming. Field observations in a subarctic tundra heath with long-term (13-year) warming treatments were extensively used for parameterizing and evaluating BVOC-related processes (photosynthesis, emission responses to temperature and vegetation composition). We propose an adjusted temperature (T ) response curve for Arctic plants with much stronger T sensitivity than the commonly used algorithms for large-scale modelling. The simulated emission responses to 2 • C warming between the adjusted and original T response curves were evaluated against the observed warming responses (WRs) at short-term scales. Moreover, the model responses to warming by 4 and 8 • C were also investigated as a sensitivity test. The model showed reasonable agreement to the observed vegetation CO 2 fluxes in the main growing season as well as day-to-day variability of isoprene and monoterpene emissions. The observed relatively high WRs were better captured by the adjusted T response curve than by the common one. During 1999-2012, the modelled annual mean isoprene and monoterpene emissions were 20 and 8 mg C m −2 yr −1 , with an increase by 55 and 57 % for 2 • C summertime warming, respectively. Warming by 4 and 8 • C for the same period further elevated isoprene emission for all years, but the impacts on monoterpene emissions levelled off during the last few years.
At hour-day scale, the WRs seem to be strongly impacted by canopy air T , while at the day-year scale, the WRs are a combined effect of plant functional type (PFT) dynamics and instantaneous BVOC responses to warming. The identified challenges in estimating Arctic BVOC emissions are (1) correct leaf T estimation, (2) PFT parameterization accounting for plant emission features as well as physiological responses to warming, and (3) representation of long-term vegetation changes in the past and the future. Peñuelas and Staudt, 2010;Possell and Loreto, 2013). BVOC synthesis is regulated by enzyme activity, and many compounds are emitted in a temperature-dependent (T ) and lightdependent (Q) manner (Li and Sharkey, 2013). BVOCs released into the atmosphere react with hydroxyl radicals (OH), which could reduce the atmospheric oxidative capacity and therefore lengthen the lifetime of methane (CH 4 ), as a potent greenhouse gas (Di Carlo et al., 2004;Peñuelas and Staudt, 2010). An increase in BVOC emissions could also elevate the tropospheric ozone (O 3 ) concentration when the ratio of BVOCs to NO x (BVOCs / NO x ) is high (Hauglustaine et al., 2005), and increase secondary organic aerosol (SOA) formation (Paasonen et al., 2013). BVOCs could also limit ozone formation when the BVOCs / NO x ratio is low, a situation in which NO x can react with O 3 (Pusede and Cohen, 2012). Global estimates of non-methane BVOC emissions are in the range of 700-1000 Tg C yr −1 , of which isoprene and monoterpenes contribute most of the emissions (∼ 70 and 11 %, respectively; Sindelarova et al., 2014). The modelled emission rates for isoprene are of similar magnitude as for CH 4 (Arneth et al., 2008). However, the current estimates of regional emission distributions are highly uncertain for both isoprene and monoterpenes for two reasons: (1) the current emission estimates are based on field studies mainly covering tropical, temperate and boreal ecosystems (Guenther et al., 2006), lacking observational data for the subarctic and Arctic, and (2) the uncertainties in driving variables (vegetation distribution and seasonality, climate and environmental data, including soil water availability and the spectrum of the incoming light, abiotic and biotic stress) and in emission responses to these drivers (Guenther et al., 2006;Arneth et al., 2008). For instance, plants adapted to the cold environment of the Arctic appear to respond to warming differently than plants from low latitudes . Until now, the emissions from high latitudes (including the Arctic and the subarctic) have been assumed to be minimal due to low foliar coverage, T and plant productivity (Guenther et al., 2006;Sindelarova et al., 2014). However, recent observations from the Arctic have indicated the need for revising the current assumption, as higher emissions from both plants and soils than anticipated in large-scale models have been measured (Ekberg et al., 2009;Holst et al., 2010;Potosnak et al., 2013;Rinnan et al., 2014;Schollert et al., 2014;Kramshøj et al., 2016). Furthermore, field experiments focusing on the effects of climate warming on BVOC emissions have found unexpectedly high responses of BVOC release to a few degrees of warming (Tiiva et al., 2008;Faubert et al., 2010;Valolahti et al., 2015;Kramshøj et al., 2016;Lindwall et al., 2016), which has underlined the potentially significant role of Arctic BVOC emissions under changing climate. The Arctic is warming at approximately twice the global rate (IPCC, 2013) and the warming-induced drastic vegetation changes (AMAP, 2012) could impose substantial changes in BVOC emissions.
Both isoprene and monoterpenes are produced through the 2-C-methyl-D-erythritol4-phosphate/1-deoxy-D-xylulose-5phosphate (MEP-DOXP) pathway and are reaction products of their chief precursors, glyceraldehyde-3-phosphate (G3P) and pyruvate. G3P is produced along the chloroplastic Calvin cycle. Mechanistic models have often linked the biosynthesis of isoprene and monoterpenes with photosynthesis processes (Niinemets et al., 1999;Martin et al., 2000;Zimmer et al., 2003;Grote et al., 2014). In the short term (hours-days), the responses to Q and T of isoprene and monoterpene production are very similar to those of photosynthesis, but with a higher T optimum for BVOC production than photosynthesis (Guenther et al., 1995;Arneth et al., 2007). Furthermore, some monoterpenes can be emitted from storage pools in plant organs, e.g. glands or resin ducts (Franceschi et al., 2005). Along with the short-term responses, the long-term (days or longer) BVOC dynamics are affected by vegetation composition changes Valolahti et al., 2015), vegetation phenology (Staudt et al., 2000;Hakola et al., 2006), past weather conditions (Ekberg et al., 2009;Guenther et al., 2012) and growing conditions, e.g. soil water and nutrient availability (Possell and Loreto, 2013), atmospheric CO 2 (Wilkinson et al., 2009) and ozone levels (Loreto et al., 2004;Calfapietra et al., 2007). Here, we use a processbased ecosystem model to represent BVOC synthesis and emissions. The model simulates vegetation composition dynamically and represents long-term growing environment effects, and thus it is useful in terms of predicting long-term emission responses to environmental changes.
Usually, estimates of BVOC responses to Q and T are based on the Guenther algorithm (referred to here as G93: Guenther et al., 1993) and observed emission rates are often standardized to emission capacity at standard conditions (T of 30 • C and photosynthetically active radiation (PAR) of 1000 µmol m −2 s −1 ) using the G93 algorithm to allow for comparison with other observations. Potosnak et al. (2013) fitted leaf-level isoprene emission rates to T and Q in moist acidic tundra and found that the G93 algorithm characterized emissions well with the T response, but not Q response. However, Ekberg et al. (2009) found that the T response of the G93 algorithm is not sensitive enough to capture the observed high T responses of wet tundra sedges, which was further supported by other studies in the high latitudes (Faubert et al., 2010;Holst et al., 2010). Furthermore, species-specific emission profiles Schollert et al., 2015;Vedel-Petersen et al., 2015) have not yet been integrated into the modelling of Arctic BVOC emissions Guenther et al., 2012;Sindelarova et al., 2014). These need to be included as a trait of plant functional types (PFTs), especially when studying the drastic impacts of climate change on vegetation composition as well as BVOC emissions in the Arctic. In addition, tundra plants with relatively dark surfaces and low growth forms (commonly less than 5 cm tall) may experience much higher leaf T than the air T at 2 m height provided by weather stations (Körner, 2003;Scherrer and Körner, 2010;Lindwall et al., 2016), which could lead to larger emissions than anticipated in current models.
The aim of this work was to integrate the observed emission features of Arctic plants into a process-based ecosystem model in order to improve the current model estimations of Arctic BVOC emissions and to advance our understanding regarding emission dynamics for Arctic ecosystems in a warming future. The process-based dynamic ecosystem model LPJ-GUESS (Lund-Potsdam-Jena General Ecosystem Simulator) (Smith et al., 2001(Smith et al., , 2014 was used as a platform to simulate short-term and long-term responses of BVOC emissions to changes in climate for Arctic plants. The model links isoprene and monoterpene production with photosynthesis (Arneth et al., 2007;Schurgers et al., 2009). For the application to a subarctic heath tundra, the process parameterization utilized field observations of long-term (13year) warming treatment effects on vegetation composition and BVOC emissions (Tiiva et al., 2008;Faubert et al., 2010;Valolahti et al., 2015). The specific objectives of this study were (1) to capture the observed T response of BVOC emissions for a subarctic ecosystem, (2) to address the importance of short-term and long-term impacts of warming on ecosystem as well as BVOC emissions, and (3) to diagnose key model developments needed to better present BVOC dynamics for the Arctic region.

Study area and observational data
The data used in this modelling study were collected at a dwarf shrub/graminoid heath tundra located in Abisko, northern Sweden (68 • 21 N, 18 • 49 E). The vegetation consists of a mixture of evergreen and deciduous dwarf shrubs, graminoids and forbs. A long-term field experiment was established at this site in 1999 to investigate the effects of climate warming and increasing litter fall, resulting from the expanding tundra vegetation, on the functioning of the ecosystem. The experiment included control (C), warming (W), litter addition (L) and combined warming and litter addition (WL) treatments . In the current study, we only focused on the observations from the C and W treatments. Each treatment, covering an area of 120 × 120 cm, was replicated in six blocks. The W treatments used opentop chambers (OTCs), which passively increased air T by around 2 • C and also caused around 10 % reduction in PAR .
During the years 2006, 2007 and 2012, BVOC emission rates were measured for all plots by sampling air from transparent polycarbonate chambers into adsorbent cartridges using a push-pull enclosure technique and analysis by gas chromatography-mass spectrometry. The enclosure covered a 20 cm × 20 cm area in each plot. The air T inside the enclosure and PAR in ambient conditions were measured during the sampling. For 2006-2007, the datasets for isoprene emission can be found in Tiiva et al. (2008) and those for monoterpenes in Faubert et al. (2010). For the year 2012, isoprene and monoterpene emissions have been published by Valolahti et al. (2015). Notably, BVOCs in this study only refers to isoprene and monoterpenes. Closed chamber-based CO 2 fluxes were measured in the same area for 2006, 2007, 2010 and 2012 (data from 2006 and 2007 were published in Tiiva et al., 2008, whilst data from 2010 and 2012 have not been published before). Species composition and coverage in the plots in the same years were estimated by pointintercept-based method, in which a hit is recorded each time a plant species is touched by a pin lowered through 100 holes covering the plot area of 20 cm × 20 cm (Tiiva et al., 2008;Valolahti et al., 2015). Species composition was measured in June for 2006, 2010 and 2012, and in June, July and August for the year 2007.

LPJ-GUESS general framework
LPJ-GUESS is a climate-driven dynamic ecosystem model with mechanistic representations of plant establishment, mortality, disturbance and growth as well as soil biogeochemical processes (Smith et al., 2001;Sitch et al., 2003). Vegetation in the model is defined and grouped by PFTs, which are based on plant phenological and physiognomic features, combined with bioclimatic limits (Sitch et al., 2003;Wolf et al., 2008). The model has been widely and successfully applied for simulating vegetation and soil carbon fluxes as well as vegetation dynamics at different spatial scales (Wolf et al., 2008;Hickler et al., 2012;Smith et al., 2014;Tang et al., 2015). In the model, individuals of each PFT in the same patch (replicate unit in the model, representative of vegetation stands with different histories of disturbance and succession) can compete for light and soil resources. Plant establishment and mortality are represented as stochastic processes, but influenced by life history, resource status and demography (Smith et al., 2014). For summergreen plants, an explicit phenological cycle is implemented, which is based on the accumulated growing degree-day (GDD) sum for leaf onset and full leaf cover.
In LPJ-GUESS, a generalized Farquhar photosynthesis model (Farquhar et al., 1980;Collatz et al., 1991) for largescale modelling is used to simulate canopy-level carbon assimilation and the generalized model is built on the assumption of optimal nitrogen (N) allocation in the vegetation canopy (Haxeltine and Prentice, 1996a, b). Daily net photosynthesis is estimated using a standard non-rectangular hyperbola formulation, which gives a gradual transition between the PAR-limited (J E ) and the rubisco-limited (J C ) rates of assimilation (Haxeltine and Prentice, 1996b). For C 3 plants, J E is a function of the canopy-absorbed PAR, the intrinsic quantum efficiency for CO 2 uptake (α c 3 ), the CO 2 compensation point ( * ) and the internal partial pressure of CO 2 (p i ) (Collatz et al., 1991;Haxeltine and Prentice, 1996b). J C is related to the maximum catalytic capacity of rubisco per unit leaf area (V m ), * , p i and the Michaelis-Menten constant for CO 2 and O 2 . Stomatal conductance influences the intercellular CO 2 , p i and canopy transpiration.

BVOC modelling
In LPJ-GUESS, isoprene (Arneth et al., 2007) and monoterpene (Schurgers et al., 2009) emissions are simulated as a function of the photosynthetic electron flux. The productions of isoprene (E I ) and monoterpenes (E M ) are computed as where J is the rate of photosynthetic electron transport and α converts photon fluxes into terpenoid units. The synthesis of both compounds is linked to J (Niinemets et al., 1999(Niinemets et al., , 2002 and a fraction (ε) of the electron transport contributing to terpenoid production (Eq. 2) is determined from a plantspecific fraction under standard conditions (ε S , usually at a T of 30 • C and a PAR of 1000 µmol m −2 s −1 ) which is adjusted for leaf T , seasonality (σ ), and atmospheric CO 2 concentration: The standard fraction ε S is computed from the often reported standard emission rate (emission capacity) together with the simultaneously estimated photosynthetic electron flux under these standard conditions (standard T and PAR) in the model. The choice of different T and PAR as standard conditions will influence the value for ε S and thus the estimated emission rate under different conditions. The T response corrects for the T optimum for terpenoid synthesis, which is higher than that for photosynthesis: The parameter α τ represents the T sensitivity and the standard temperature (T S ) is often 30 • C (adjusted to 20 • C in this study). In the model, daily mean T (T d , model input) has been adjusted to daylight hours T based on day length as well as daily T range (Arneth et al., 2007) and the daytime T is used for calculating daily emission rates. For the study in the subarctic, the often-used reference T S of 30 • C and the T responses (α τ ) were adjusted based on the observation data and will be discussed below. The seasonality function, f (σ ), was applied to both isoprene and monoterpene production and is based on a degree-day sum in spring and daylength thresholds in autumn (Arneth et al., 2007;Schurgers et al., 2009). The atmospheric CO 2 concentration enhances terpenoid synthesis when the concentration is lower than ambient, and vice versa, which is represented by the function f (CO 2 ) (Arneth et al., 2007). The model assumes that both isoprene and monoterpenes are produced in the same pathway and that they respond to CO 2 concentration in the same way. For monoterpenes, a storage pool (m) is assigned to represent the specific (long-term) storage of monoterpenes within a leaf (Schurgers et al., 2009). The storage pool is only implemented for coniferous and herbaceous PFTs (see Table S1 in the Supplement). The emission of monoterpenes from the storage (E Ms ) is a function of T d and m with an average residence time (τ ). τ S is the residence time at the standard T of 30 • C (adjusted to 20 • C in this study, consistent with the modification on the T responses of terpenoid synthesis). The residence time τ is adjusted based on the standard condition τ S for T d responses with a Q 10 relationship.
In LPJ-GUESS, the BVOC response to light resides in the photosynthesis processes (light dependence of J in Eq. 1). Additionally, considering the high sensitivity of BVOC production to leaf T , the model applies a computation of leaf T based on air T and energy balance constraints (Arneth et al., 2007;Schurgers et al., 2009). The calculation of leaf T in the model was based on solving the leaf energy balance, where the incoming shortwave and longwave radiation are balanced by the outgoing longwave radiation and sensible heat fluxes as well as latent heat loss. The existing leaf energy balance equations appeared to underestimate the incoming longwave radiation under overcast conditions, which has been updated by specifically considering the cloud emission of longwave radiation relative to clear-sky condition (Sedlar and Hock, 2009). The estimated leaf T , rather than air T , was used for both photosynthesis and BVOC synthesis. Water loss (latent heat fluxes) is regulated by stomatal conductance and soil water content, which is also linked to leaf T estimation in the model.

Simulation setup 2.3.1 Input data
The daily climate data of air T , air T range and precipitation for the period 1984-2012 (Callaghan et al., 2013;Tang et al., 2014) were provided by the Abisko scientific research station (Abisko Naturvetenskapliga Station, ANS). Four gaps in daily radiation data from ANS (during the periods of 1 January-30 June 1984, 9-16 June 2016, 13-15 February 2007 July-17 August 2011) were filled with the Princeton reanalysis dataset (Sheffield et al., 2006) for the grid cell nearest Abisko. The annual CO 2 concentrations for the whole study period  were obtained from McGuire et al. (2001) and TRENDS (http://cdiac.esd.ornl.gov/). The air T inside the enclosure and ambient PAR at canopy level were also used as the model inputs for each measuring day (Tiiva et al., 2008;Faubert et al., 2010;Valolahti et al., 2015).

Plant functional types
The dominant plant species from the observations  were divided into seven PFTs (Table 1). The PFT parameters (see Table S1) were mainly derived from previous studies for the Arctic region using LPJ-GUESS (Wolf et al., 2008;Miller and Smith, 2012;Tang et al., 2015), but the Arctic PFT lists were extended to consider BVOC emission characteristics. The low summergreen shrubs (LSS) were divided into a Salix type (SLSS; high isoprene emitter) and a non-Salix type (NSLSS; e.g. Betula nana dominance, predominantly monoterpenes rather than isoprene emitters) (Schollert et al., 2014;Vedel-Petersen et al., 2015). Furthermore, due to the abundance of prostrate dwarf shrubs (PDS) in the study area, distinguishing PDS (canopy height lower than 20 cm) from low shrubs (canopy height lower than 50 cm) was implemented through adjusting parameters controlling vegetation height. The PDS type was further divided into two PFTs with evergreen and deciduous phenology. Moss, widely appearing in the study area, was not dis-tinguished from forbs and lichens, due to limited data for parameterizing moss physiognomic features and their preferable growing conditions. In LPJ-GUESS, the crown of each tree is divided into thin layers (original value is 1.0 m in a forest canopy) in order to integrate PAR received by each tree. The thickness of this layer was reduced to 10 cm in this study to better capture the vertical profile of low and prostrate shrubs. In addition, the original specific leaf area (SLA, m 2 kg C −1 ) values in LPJ-GUESS were estimated based on a fixed dependency on leaf longevity (Reich et al., 1997). In our study, a fixed SLA was assigned to each PFT (Oberbauer and Oechel, 1989) to improve the simulated leaf area index (LAI) for Arctic plants. Emission capacities for the PFTs were determined from available leaf-level measurement data from the subarctic and Arctic. The details about the data sources for parameterizing emission capacity at 30 • C (E IS30 , E MS30 ) and 20 • C (E IS20 , E MS20 ) can be found in Table S2 and the averaged emission capacities (among all literature data in Table S2) for each PFT as well as the representative plant species can be found in Table 1. The emission rates from the literature are generally provided as standardized emission capacities at 30 • C using the G93 algorithm, and these values were further Figure 1. The observed isoprene emission rates in relation to the chamber air temperature in July over three field seasons (2006,2007,2012) in the Abisko tundra heath. rescaled to 20 • C using the adjusted T response curve from this study (Fig. 1).

Model calibration and evaluation
The modelled CO 2 fluxes, LAI and BVOC T response were first calibrated before evaluating the modelled daily BVOC emission rates. Two out of four years' (2006 and 2007) measured net ecosystem production (NEP), ecosystem respiration (ER) and estimated gross primary production (GPP) as well as point-intercept-based species composition were used for calibrating. The data for the other 2 years (2010 and 2012) were used for evaluating the simulated carbon cycle processes. Previous studies focusing on light responses of NEP for Arctic plants (Shaver et al., 2013;Mbufong et al., 2014) have reported relatively low quantum efficiencies (α c 3 ) caused by overall low sun angle conditions and low leaf area. A thorough sensitivity study of parameters used in LPJ-GUESS (Pappas et al., 2013) has found that α c 3 is the most influential parameter in terms of the simulated vegetation carbon fluxes. Also, a pre-evaluation of the modelled CO 2 fluxes with the observations in this study using the default α c 3 value (0.08) has found a large overestimation of both GPP and ER (not shown). Therefore, a sampling of α c 3 (using the range of 0.02 to 0.125 µmol CO 2 µmol photons −1 , proposed by Pappas et al., 2013) was conducted to find the best value to depict the observed GPP, ER and LAI of the years 2006 and 2007 for the subarctic ecosystem (Fig. S1 in the Supplement). After calibration, the model was evaluated with the simulated CO 2 fluxes and vegetation composition using the observed CO 2 fluxes and the point-intercept-based plant coverage data from 2010 and 2012, respectively.
The daytime air T in the study area is often below 20 • C (Ekberg et al., 2009), and standardization of terpenoid emissions to 20 • C, instead of 30 • C, has been suggested for modelling in boreal and Arctic ecosystems (Holst et al., 2011, Ekberg et al., 2009 due to plant adaptation to low T environment. In the model, the photosynthetic electron fluxes under standardized conditions are simulated in order to convert the input emission capacity to the standard fraction (ε S , see Eq. 2). The choice of the standardized T (used in Eq. 3 as well as in estimating photosynthesis rates at this T ) will influence the estimated fraction of electron fluxes for BVOC synthesis. In this study, data fitting to the suggested standard T of 20 • C was conducted using the observed ecosystemlevel isoprene emission rates in July together with measurement chamber air T from the C plots. The observations were mostly conducted during daytime with relatively high PAR values, and therefore the response of the emission rates to light was not specifically considered in the current data fitting. Potential feedbacks from the variations in the atmospheric CO 2 concentration were ignored for the 3 years with isoprene sampling (a rough model estimation of ∼ 3% reduction in emissions between 2006 and 2012). The data collected from different blocks were separated for the curve fitting and the parameters controlling T response (α τ in Eq. 3) were determined (Fig. 1). An adjusted α τ value of 0.23 was chosen after fitting all the data from July over 3 years' measurements. Apart from the low R 2 value for block 1, the data were well captured by the exponential shape (R 2 ≥ 0.8) of the T response curve. The calibrated T responses were used for standardizing leaf-level emission rates (see M IS20 and E MS20 , Table 1) as well as estimating emission rates in the model. This adjusted T response was also evaluated with the observed enclosure air T and monoterpene emission rates in July (R 2 = 0.66 for all blocks). The abundance of each PFT was evaluated using simulated LAI against the point-intercept-based vegetation composition. The species were grouped into the corresponding PFTs for comparison and the point-intercept-based hits within the same PFT group were summed. The summed hits were divided by 100 pin hits to compare with the modelled LAI. The point-intercept-based species abundances and LAI are not comparable one to one throughout growing seasons, since the measurement could include pin hits on different plant parts, whereas LAI only explains leaf coverage. However, the point-intercept-based coverage approaches leaf coverage when the deciduous leaves become fully developed during the growing season.
After calibration of the modelled CO 2 fluxes and LAI, the modelled isoprene and monoterpene emission rates were compared with the observations. The simulated daytime average emissions (µg C m −2 h −1 , daytime emission rates divided by day length) do not allow an accurate comparison with the observed emission rates, which were typically obtained in the middle of the day (between 09:00 and 17:00). Therefore, an additional estimate of the emission rates for the conditions prevailing during the sampling was made. This was done by computing the emission applying the measured air T inside the enclosure and PAR during the sampling time for photosynthesis and BVOC emissions. This computation was performed twice: once using the original T response (α τ = 0.1, T S = 30 • C, E IS30 and E MS30 , Eq. 3) and once with the adjusted T response (α τ = 0.23, T S = 20 • C, E IS20 and E MS20 , Eq. 3 and Fig. 1).
The model's performance in modelling BVOC emissions was evaluated by Willmott's index of agreement (A) (Eq. 5) and mean bias error (B) (Eq. 6). The index A describes the agreement between the modelled fluxes (E i ) with the observed (O i ) and a value close to 1 indicates a good agreement. The index B estimates the mean deviation between the modelled and observed values (Willmott et al., 1985) and values close to 0 indicates models' good agreement with observations.
where O is the observed mean value and N is total number of data records.

Effect of warming
To simulate the observed warming responses from the OTCs, a warming of 2 • C was imposed in the model for the growing season (the period with OTC warming) (Tiiva et al., 2008;Valolahti et al., 2015). The modelled warming responses (WRs, difference between C and W treatments) using the original T response and the adjusted T response were compared with the observed WRs. Furthermore, additional simulations with a warming by 4 and 8 • C, reflecting the range of climatic projections in this region (IPCC, 2013), were also conducted to test for the anticipated ecosystem-scale responses to different levels of warming.

Modelled CO 2 fluxes and vegetation composition
The simulated ecosystem CO 2 fluxes and LAI were sensitive to the parameter value chosen for α c 3 , which describes the efficiency in converting solar radiation to carbohydrates, and which was varied between 0.02 and 0.125 µmol CO 2 µmol photons −1 following Pappas et al. (2013) (Fig. S1). For CO 2 fluxes, the lowest rootmean-square error (RMSE) values occurred at 0.035 µmol CO 2 µmol photons −1 for GPP and ER, while the lowest RMSE value for LAI was 0.051 µmol CO 2 µmol photons −1 when comparing with the observations for 2006 and 2007. A value of 0.040, consistent with the study by Shaver et al. (2013), was selected for α c 3 to limit the RMSE values of the modelled CO 2 fluxes and LAI. Using this value for α c 3 , the model captured the observed day-to-day variations as well as the magnitude of the chamber-based GPP, ER and NEP for 2010 and 2012, with an overestimation of CO 2 fluxes (particularly for the early growing seasons, Fig. 2) and a large underestimation of LAI (Fig. 3). For the year 2012, the model showed large overestimations of the observed GPP and ER for the limited number of measurements in this growing season. For the five PFT groups, the modelled growing season LAI values for 2010 and 2012 were much lower than the pointintercept-based coverage estimations from the field observations (note different left and right axis scales in Fig. 3 to allow comparison of relative changes in response to warming), except for the Salix-type summergreen shrubs and deciduous prostrate dwarf shrubs (SLSS + SPDS). The dominance of two vegetation groups in the C plots -forbs/lichens and evergreen shrubs -was consistent between the modelled and the observed.
In response to 2 • C warming, the modelled LAI for the shrub PFTs (SLSS + SPDS, NSLSS, LSE + EPDS) showed an increase, while the modelled LAI for graminoids and forbs/lichens largely decreased (Fig. 3). For the two groups of shrubs (NSLSS and LSE + EPDS), the modelled increase is in agreement with the observations. However, the observed large increase in the coverage of forbs/lichens and a decreased coverage of graminoids in the W treatments for the year 2010 and 2012 were not captured by the model.

Modelled BVOC emissions
BVOC emissions are closely linked to leaf as well as ecosystem development. Simulating seasonal variation in leaf area and vegetation composition enables us to assess the model performance in representing short-term emission changes in response to T and PAR, as well as long-term changes in vegetation development and distribution. The seasonal variations in the modelled daily BVOC emissions as well as the span of all BVOC samplings over 3 years are presented in Fig. S2.

Daily emissions Emission rates in the control (ambient) conditions
The observed air T and PAR showed day-to-day variations through the sampling periods (Fig. 4e), which resulted in strong daily variations in the observed BVOC emissions ( Fig. 4a and c). These observed variations in isoprene and monoterpene emissions were generally captured by the model for 2006 and 2007. For the year 2012, the model overestimated both isoprene and monoterpene emission rates over the three sampling days. Noticeably, the model used air T at 2 m height from the ANS station to extrapolate the leaf T for estimating daily BVOC emissions (Fig. S2), while the observed air T and PAR during the sampling hours were used for modelling the emissions to directly compare with the observed (Fig. 4). The modelled high emission rates for a few days (e.g. 10 July 2007, 14 June 2012) were directly linked to the observed high T and PAR at the canopy level (Fig. 4e). Averaged over all measuring days in 2006 and 2007, the modelled and observed isoprene emission rates were 46.6 and 34.7 µg C m −2 h −1 , and the modelled and observed monoterpene emission rates were 8.5 and 5.3 µg C m −2 h −1 , respectively. For the year 2012, the modelled emission rates (80.4 and 14.9 µg C m −2 h −1 for isoprene and monoterpenes, respectively) were much higher than the observed (9.1 and 0.5 µg C m −2 h −1 , for isoprene and monoterpenes, respectively). The large overestimation by the model in the year 2012 was also seen for GPP and ER (Fig. 2).

Emission responses to 2 • C warming
In response to warming by the OTCs, the observed enclosure air T in the W plots was 2.1 • C higher than that in the C plots averaged over the three growing seasons with observations. For isoprene, the observed magnitudes of WRs (Fig. 4b) were captured reasonably well by the model, ex-  cept for 5 August 2007. For this day, the air T in the W was higher than in the C plots, but the PAR value was lower in the W than in the C plots (Fig. 4e). Averaged over 3 years, the simulated and observed isoprene WRs were 19.6 and 28.4 µg C m −2 h −1 , respectively. Warming increased the observed isoprene emissions by 95 % but only increased the modelled emissions by 37 % (dividing the averaged WRs by the averaged emissions for the days on which measurements were made). For monoterpenes, the modelled and observed WRs were 6.1 and 4.0 µg C m −2 h −1 , respectively. Averaged over three growing seasons, warming increased the observed monoterpene emissions by 93 %, and the modelled emission by 63 % (dividing the averaged WRs by the averaged emissions for the days on which measurements were made).
These modelled WRs obtained with the adjusted BVOC T response (α τ = 0.23, T S = 20 • C, Eq. 3) were further compared with the simulation using the original T response (α τ = 0.1, T S = 30 • C, Eq. 3). For isoprene (Fig. 5a), the simulation using the adjusted T response showed a substantial increase in the modelled WRs as well as a better agreement with the observations (A = 1.16, B = −8.85) than the simulation using the original T response (A = 1.47, B = −27.26). The modelled WRs using the original T response largely underestimated the observed high WRs. Averaged over 3 years, the isoprene WRs modelled using the original T response (used at a global scale) only gave 4 % of the observed WRs, while the WRs modelled using the new T response captured 69 % of the observed WRs (using the modelled average WR to divide with the observed average WRs). For monoterpenes, the WRs modelled using the adjusted T response (A = 0.80 and B = 2.13) showed a moderate improvement as compared to using the original T response (A = 1.35 and B = −2.83). The modelled WRs using the original T response underestimated the observed WRs by 72 %, but the modelled WRs using the adjusted T response overestimated the observed WRs by 53 %. For the year 2007, the observed high monoterpene WRs was better captured by the simulated WRs with the adjusted T response. As for the modelled emission rates, the overestimation of the observed WRs also mainly occurred in 2012.

Annual emissions
A comparison of the simulated annual BVOC emissions from the C and W treatments demonstrated that the 2 • C warming during the growing seasons increased both isoprene and monoterpene annual emissions. Averaged over 13 years, this warming increased annual isoprene and monoterpene emissions by 55 and 57 %, respectively (p<0.01, Mann-Whitney test). The modelled emissions showed strong interannual variations in response to warming (Fig. 6). For the warmest year (2011), the W treatment increased annual isoprene and monoterpene emissions by 99 and 94 %, respectively. The mean annual isoprene and monoterpene emissions in the C for 1999-2012 were 20 and 8 mg C m −2 yr −1 , respectively. For the 3 years with BVOC sampling, the modelled average WRs were 58 and 70 % for annual isoprene and monoterpene emissions, respectively. The modelled annual WRs were of similar magnitude as the modelled daily average WRs (data not shown) for the days with BVOC samplings (63 % for isoprene and 81 % for monoterpenes).
The simulations imposing the warming by 4 or 8 • C during the same period as the 2 • C warming increased annual isoprene emissions by 120 and 247 %, respectively (p<0.01, Mann-Whitney test) and annual monoterpene emissions by 87 and 167 %, respectively (p<0.01, Mann-Whitney test). For isoprene, the strongest WRs of all levels of warming appeared in 2011. Higher levels of warming further elevated isoprene emissions for all years, but the impact on monoterpene emissions levelled off due to a decreasing coverage of evergreen prostrate dwarf shrubs (EPDS) with 8 • C warming. The decrease in coverage of EPDS only occurred for the last few years with 4 • C warming. The different levels of warming generally increased shrub growth, but largely decreased the coverage of forbs/lichens and graminoids (CLM and GRT) (data not presented). At annual scale, the long-term vegetation changes associated with warming by 4 or 8 • C showed strong impacts on BVOC emissions.

Emission rates
The modelled day-to-day variations in ecosystem CO 2 fluxes (Fig. 2) and BVOC emissions (Fig. 4) generally followed the observations, in spite of the poor representation of the observed vegetation composition (Fig. 3). The mismatch between the modelled LAI and the observed vegetation coverage is likely partly due to the fact that LAI only includes the areal coverage by leaves, whereas the point-intercept-based vegetation coverage also includes coverage detected of other aboveground plant parts, like stems. Further, the mismatch may also be caused by an underestimation of the allocation of assimilated carbon to foliage in LPJ-GUESS and/or too low SLA values (Table S1). In LPJ-GUESS, the carbon allocation among different living tissues follows four allometric equations to control the structural development of each modelled plant individual (see Eqs. 1-4 in Sitch et al., 2003). The allometric parameters for some of the Arctic PFTs used in this study were validated by Wolf et al. (2008) derived for a model applying a quantum efficiency α c 3 of 0.08 at the regional scale, which may require further justification after the reduction in α c 3 that was applied here to match the observed daily CO 2 fluxes. The reduced quantum efficiencies reflect the growth environment with low T and low sun angle in high latitudes (Shaver et al., 2013), but more observations are still needed to better quantify light use efficiency of Arctic plants (Dietze et al., 2014). Furthermore, Van Wijk et al. (2005) found a close linkage between total foliar N content and LAI for Arctic plants, which was further supported by Campioli et al. (2009) for an Arctic ecosystem dominated by Cassiope tetragona. However, the current simulations neither include C-N interactions nor consider potential impacts of N limitation on plant development (Smith et al., 2014), which need to be improved in future model simulations in this region . The subdivision of Arctic PFTs into smaller groups to specifically consider isoprene and monoterpene emission features was shown to be important for capturing the emission dynamics in this heath tundra ecosystem. The development of parameterizations for Arctic PFTs also requires considering the phenological and physiognomic features of mosses (currently aggregated in the CLM-type PFT, Table S1), which may bring additional uncertainties to the modelled LAI. The current evaluation of the modelled LAI with the point-intercept-based measurements of plant coverage cannot disregard uncertainties from the field method itself, such as subjective judgement of species from each hit, and sampling inclining angles (Wilson, 2011). Also, the seasonal variation in leaf development as well as the randomly selected blocks from the hetero-geneous landscape may further complicate the comparison of the simulated LAI with the local observations. Capturing the start of the growing season in the model is also crucial for depicting the dynamics of seasonal CO 2 fluxes (Tang et al., 2015). The overestimated GPP at the beginning of growing seasons (Fig. 2a) suggests uncertainties in modelling the time of its start. The current algorithm for detecting the start of the growing season in large-scale applications (Sykes et al., 1996) may not be sensitive enough for prediction of budburst of Arctic plants (Pop et al., 2000).
The modelled annual isoprene and monoterpene emissions, 20 and 8 mg C m −2 yr −1 for 1999-2012, correspond to less than 0.1 % of the modelled GPP. The modelled emission rates are not only linked to the modelled photosynthesis fluxes but also determined by the emission capacity assigned to each PFT (see Tables 1 and S2). For some PFTs (e.g. the Salix-type and prostrate summergreen shrubs, SLSS and SPDS), the emission capacities in Table 1 are of similar magnitude to observed values that are applied in large-scale models for boreal forests (see Table 2 in Rinne et al., 2009). The observed relatively low emissions in comparison with lower latitudes Sindelarova et al., 2014) are mainly caused by low T and plant biomass, and not by low emission capacities (Holst et al., 2010).
The numbers for the estimated annual emissions are still highly uncertain, considering the dissimilarities to the observations in the modelled LAI, early season CO 2 fluxes and the overestimation of daily isoprene and monoterpene emissions of a few days. The observed low values of CO 2 fluxes (GPP and ER) and BVOC emissions in 2012 could be due to harmful effects of an insect outbreak in the nearby birch forest (Hanna Valolahti, personal observation). The potential impacts from insect outbreaks have not been explicitly included in the model. When both T and PAR were high (e.g. on 6 July 2007), the model tended to overestimate the emission rates, which could suggest that the stronger T sensitivity that was obtained in this study does not extend to these high temperature values. Furthermore, the estimated emission rates may be more robust for isoprene than for monoterpenes, because (1) the adjusted T response curve was only applied for monoterpene production, and there is a lack of data for evaluating T responses of monoterpene emissions from storage pools (Eq. 4), and (2) there are more studies supporting CO 2 inhibition on isoprene emissions (Arneth et al., 2007) than on monoterpenes (Peñuelas and Staudt, 2010). Therefore, more laboratory experiments in controlled conditions testing BVOC responses (especially monoterpenes) of Arctic plants to different environmental variables could largely reduce the abovementioned uncertainties. Based on the current estimation, the relative magnitudes of isoprene and monoterpene emissions from this site may not contribute significantly to the global number. However, the highly reactive compounds emitted by plants could undergo chemical reactions in the local/regional atmosphere and provide feedbacks to the climate. Furthermore, the warming-induced strong increase in emissions could indicate an increasing role of BVOCs in the local atmospheric chemistry and also global emission magnitudes for future conditions.
Relative to isoprene emission, the magnitude of monoterpene emissions was much lower since the species in the study area were mostly considered to be isoprene emitters (Tiiva et al., 2008;Faubert et al., 2010). The observed monoterpene emissions were generally low for the sampling days (see Fig. S2), which could bias the evaluation. More observations in the higher T range would enhance our confidence in the new T response function, specifically for monoterpenes. Furthermore, the current observations of BVOC emissions only covered the main growing season. Sampling over a longer season would help to improve the parameterization of the partitioning over direct emission and storage, as well as the T response of emission rates from storage pools. Furthermore, ongoing 13 C labelling experiment focusing on Arctic mesocosms (Ghirardo et al., unpublished) could also help to identify the fraction of monoterpene emissions from production or storage.
The push-pull enclosure technique used for BVOC emission measurements can bring uncertainties to the measurement data: the choice of sampling time and flow rates influences temperature and humidity inside the enclosure and this, in addition to potential gas concentration changes within the enclosure, may impact the plant physiological status. The impacts also depend on the ecosystem emission rate (Niinemets et al., 2011) and sampling time of a day, considering the strong diurnal dynamics of BVOC emissions in the Arctic (Lindwall et al., 2015). The model evaluation using these half-hour-long samplings cannot avoid the influence of changed conditions inside the enclosure and of plant adaption to these conditions.

Responses to warming
The modelled increase in shrub coverage in response to the W treatment mostly followed the observations  and is consistent with the general trend in the Arctic (Wahren et al., 2005;Elmendorf et al., 2012). However, the observed increase in bryophytes is rather site-specific and was not captured by the model. In contrast, the modelled Winduced decreased coverage of graminoids and forbs/lichens agrees well with the large-scale trend identified by Elmendorf et al. (2012), who conducted a global synthesis of 61 tundra warming experiments. The decreasing soil moisture in W treatments (excluding wet ecosystems) is one of the main constraints on bryophyte coverage .
Along with vegetation community alterations, the shortterm T responses of the vegetation are central for accurately depicting daily BVOC emission responses to the W treatment. Through adjusting the BVOC T sensitivity (from α τ = 0.1, T S = 30 • C to α τ = 0.23, T S = 20 • C in Fig. 1), the simulated BVOC WRs (19.6 µg C m −2 h −1 for isoprene and 6.1 µg C m −2 h −1 for monoterpenes) became comparable to the observed responses (28.4 µg C m −2 h −1 for isoprene and 4.0 µg C m −2 h −1 for monoterpenes). The adjusted T response curve represents subarctic plants' isoprene emission responses to warming better than the original curve which has been parameterized for global simulations (Fig. 5). It further supports the earlier suggested stronger T sensitivity of BVOC emissions from Arctic plants compared to plants from other regions (Ekberg et al., 2009;Holst et al., 2010;Rinnan et al., 2014;Kramshøj et al., 2016). The commonly used T response in Guenther's algorithm (Guenther et al., 1993) is based on the Arrhenius-type dependence of enzyme activities with an optimum T around 40 • C, and the shape of the Guenther's response is very close to the exponential curve with an α τ value of 0.13 (using standard T of 30 • C) when leaf T is lower than 30 • . The high α τ value found in this study indicates that a slight T increase during summertime could cause a large increase in isoprene and monoterpene emissions from the studied cold subarctic ecosystem (Faubert et al., 2010;Holst et al., 2010). Furthermore, the adjusted T response is based on the data fitting of the observed canopy air T with hourly isoprene emission rates, and this response is used to estimate both the emission rates at sampling hour and daytime emissions in the model. The different temporal resolution for estimating daytime emissions calls for further adjustment of this T response for Arctic plants.
The underestimation of strong isoprene WRs on 5 August 2007 (157.8 µg C m −2 h −1 ) cannot be directly linked to the T and PAR differences between the C and W plots during the sampling time. The modelled emission at the C plot was 24 % lower than the observed, caused by slightly different meteorological conditions during the sampling, but the modelled WRs was 74 % lower than the observed on this date. The observed strong WRs could be linked to strong elevation of leaf T . The low-statured plants in dry to mesic tundra ecosystems are efficient in absorbing heat and thus prone to have a high leaf T on a sunny day (Schollert et al., 2014;Lindwall et al., unpublished). This can directly elevate BVOC emissions and WRs (Lindwall et al., 2016), and decouples leaf T from 2 m air T (Körner, 2003;Lindwall et al., 2016). Furthermore, for regions with underlying permafrost (not the case in this study site) in the Arctic, the potentially low ecosystem evapotranspiration can increase both ground and leaf T . Also, plants acclimated to a cold environment may drive larger emission responses once they are exposed to warmer T . The observed strong WRs can also be partly due to the potential side effects of the OTCs in the W treatment, e.g. reduced wind speed (De Boeck et al., 2012), drying of soil surfaces and increased frequency of high-temperature events . At annual to decadal timescales, the warming in the experimental plots caused changes in total plant biomass and species coverage which were found to contribute to the increase in BVOC emissions after 13 years of treatments . These indirect effects on BVOC emissions were not yet identified after 7-8 years of warming in 2006 and 2007 (Tiiva et al., 2008;Faubert et al., 2010), which highlights the importance of accurately representing the temporal dynamics of vegetation as a driver of BVOC emissions. The modelled annual emissions in response to different degrees of warming (Fig. 6) clearly elucidated the combined effects of the direct responses to summer warming with the indirect responses from vegetation changes, although the model still has limitations in representing the observed vegetation composition in detail (Fig. 3). Furthermore, these combined effects also suggest a non-linear response of BVOC emissions to different levels of warming.

Suggestions for further work
For extrapolating the current model developments to largescale (regional) applications, we suggest addressing the following issues: 1. The emission responses to T of Arctic plants could be further tested based on laboratory experiments in controlled conditions.
2. The strong decoupling of leaf T from air T and the strong dependence of BVOC emissions on leaf T (Lindwall et al., 2016) point to a need for accurately capturing leaf T in models. Long-term parallel observations of both leaf and air T will be useful for the algorithm development focusing on Arctic vegetation .
3. The subdivision of the existing PFTs into groups featuring isoprene and monoterpene emissions is encouraged for other relevant modelling studies (Grote et al., 2014), and additional data may be required for characterizing the new subgroups, such as bioclimatic limitations.
4. The potential impacts of seasonal dynamics of vegetation as well as phenology on emission capacities should be further identified with whole-season BVOC sampling (Staudt et al., 2000).
5. The responses and/or acclimation of Arctic PFTs to warmer climate should be better parameterized in the model to improve the representation of long-term vegetation effects on BVOC emissions.

Conclusions
This study has demonstrated the model's ability to depict the observed isoprene and monoterpene emission rates as well as daily variations in the BVOC emission of a subarctic tundra ecosystem. The modelled warming responses using a curve adjusted for a stronger T response showed good agreements with the observations, especially for the days with the observed strong emission responses to warming. Short-term underestimation of the observed peak of WRs was most likely linked to the underestimated leaf T during the daytime. In the long-term (days-years), a mismatch in the modelled vegetation composition could also bring uncertainty in the simulation of emission responses to warming. The model estimated the mean annual isoprene and monoterpene emissions to be 20 and 8 mg C m −2 yr −1 , with an around 55 and 57 % increase in annual emissions in response to a 2 • C warming for the period 1999-2012. For the warmest year, the 2 • C warming during the growing season resulted in 99 and 94 % increase in isoprene and monoterpene emissions. These strong warming responses of Arctic BVOC emissions have hitherto not been specifically described in large-scale models and are therefore suggested to be included, especially in estimating regional emissions from the pan-Arctic.
6 Data availability The model-simulated output used in this manuscript is available at Pangaea (doi:10.1594/PANGAEA.869465).
The Supplement related to this article is available online at doi:10.5194/bg-13-6651-2016-supplement.