Positive feedback mechanism between biogenic volatile organic compounds and the methane lifetime in future climates

A multitude of biogeochemical feedback mechanisms govern the climate sensitivity of Earth in response to radiation balance perturbations. One feedback mechanism, which remained missing from most current Earth System Models applied to predict future climate change in IPCC AR6, is the impact of higher temperatures on the emissions of biogenic volatile organic compounds (BVOCs), and their subsequent effects on the hydroxyl radical (OH) concentrations. OH, in turn, is the main sink term for many gaseous compounds including methane, which is the second most important human-influenced greenhouse gas in terms of climate forcing. In this study, we investigate the impact of this feedback mechanism by applying two models, a one-dimensional chemistry-transport model, and a global chemistry-transport model. The results indicate that in a 6 K temperature increase scenario, the BVOC-OH-CH4 feedback increases the lifetime of methane by 11.4% locally over the boreal region when the temperature rise only affects chemical reaction rates, and not both, chemistry and BVOC emissions. This would lead to a local increase in radiative forcing through methane (ΔRFCH4) of approximately 0.013 Wm−2 per year, which is 2.1% of the current ΔRFCH4. In the whole Northern hemisphere, we predict an increase in the concentration of methane by 0.024% per year comparing simulations with temperature increase only in the chemistry or temperature increase in chemistry and BVOC emissions. This equals approximately 7% of the annual growth rate of methane during the years 2008–2017 (6.6 ± 0.3 ppb yr−1) and leads to an ΔRFCH4 of 1.9 mWm−2 per year.


INTRODUCTION
Methane is the second most abundant and important noncondensable greenhouse gas with an atmospheric growth rate of +6.6 ± 0.3 ppb yr −1 (∼0.35% yr −1 ) during the years 2008-2017 [1][2][3] . Although global anthropogenic emissions of methane represent only 3% of the global CO 2 anthropogenic emissions in units of carbon mass flux, the increase in atmospheric methane concentrations has contributed to ∼23% (∼0.62 Wm −2 ) of the additional radiative forcing accumulated in the lower atmosphere since 1750 4 . The reasons for the intense growth of methane during the last decade are unclear and relate mainly to our lack of knowledge associated with globally varying emission sources (biogenic, thermogenic and pyrogenic) 5 , and the high uncertainties in quantifying the concentrations of the main sink term of methane: the hydroxyl radical 6 . According to a recent study, based on a multi-model ensemble 7 , the OH concentration variability needs to be accurately captured in order to avoid erroneous conclusions about the causes of the observed CH 4 concentration changes.
Boreal forests cover approximately one-third of all the world's forested land (~12 × 10 6 km 2 ) 8,9 . In these forests and other rural areas, one important sink term for the hydroxyl radical is the reaction with biogenic volatile organic compounds emitted from trees, understory, and soil 10,11 . The main BVOCs are the isoprenoids (isoprene, monoterpenes and sesquiterpenes), which have atmospheric lifetimes from minutes to days 12,13 .
At the Station for Measuring Ecosystem-Atmosphere Relations (SMEAR II, Hyytiälä, Finland, https://www.atm.helsinki.fi/SMEAR/ index.php), extensive research has been performed on the ecosystem-atmosphere exchange processes of BVOCs, showing that monoterpenes are the BVOCs with the highest concentrations 14,15 and OH-reactivity 16 , while isoprene only contributes 10-20% to the total VOC mass and OH-reactivity. This contrasts remarkably with the overall global emission inventory, where isoprene is by far the dominant BVOC. The emissions of isoprenoids vary depending on biological (e.g., plant species, plant-specific emission capacity, phenology, biotic and abiotic stresses) and physical (e.g., temperature, light and water availability) 17 factors. In all species the direct temperature response of BVOC emissions is exponential, but the coefficient (beta) varies with season, species and between different BVOC molecules [18][19][20] .
In this study, we investigated how projected temperature increases will affect the atmospheric concentrations of BVOCs in the Northern hemisphere, with a special focus on central-southern Finland. Furthermore, we studied how this will impact on the hydroxyl radical concentrations, and lifetime and concentrations of methane. We applied the one-dimensional chemistry-transport model SOSAA (model to Simulate the concentrations of Organic vapours, Sulphuric Acid and Aerosols) 16,21,22 and the global chemistry transport model TM5 23 to investigate the impact of higher temperatures on the emission rates and the chemistry of BVOCs.

SOSAA and TM5 simulations
With SOSAA we performed 11 simulations for the year 2016, including the control run denoted BASE. Three scenarios were performed with the temperature in (only) the emission module increased by 2, 4 and 6 K, denoted E2C0, E4C0 and E6C0, respectively. Three scenarios were performed with the temperature in (only) the chemistry module increased by 2, 4 and 6 K, denoted E0C2, E0C4 and E0C6, respectively. Three scenarios were performed with the temperature in both the emission and the chemistry modules increased by 2, 4 and 6 K, denoted E2C2, E4C4 and E6C6, respectively. Finally, one run (E6CS) was performed where the temperature was increased by 6 K in the emission module, and in the initial reactions of the BVOCs with the main oxidants OH, O 3 and NO 3 , but not in the rest of the chemistry module. The E6CS scenario was performed to investigate how the increased temperature will affect the OH concentration through OH-recycling by comparison to E6C6. More background information on the expected temperature dependencies in the chemical reaction scheme and OH-recycling mechanism are provided in the Supplementary Discussion. To ensure that meteorological processes do not affect our results, we only increased the temperatures in the emission and/or chemistry modules-all other parameters were left unchanged.
Additionally, we applied the calculated emission rate increase of monoterpenes from SOSAA in the E6C0 scenario in the global chemistry transport model TM5 for the whole boreal region to study the BVOC-OH-CH 4 feedback loop for the Northern hemisphere. By choosing this setup, we are able to quantify the direct impact of a temperature dependent BVOC concentration increase on the main methane sink term, the hydroxyl radical, in different future climate scenarios. Figure 1a, b show the concentrations of monoterpenes for the BASE run of 2016, and the differences between E6C6 and BASE at SMEAR II calculated by SOSAA. The predicted increase of monoterpene concentrations in our case scenarios originates from a rise in emission rates of 14%, 31% and 52% for the E2C2, E4C4 and E6C6 scenarios relative to the BASE run, respectively. This leads to an average increase in monoterpene concentrations of 16% if the temperature was increased by 2 K, rising to 36% and 60% for the E4C4 and E6C6 simulations, respectively. The concentrations are averaged from the surface to 3 km and the increase of monoterpene concentrations for all case scenarios are relative to the BASE run. The monoterpene concentrations of the BASE simulations averaged for daily concentrations for the four seasons in year 2016 are provided in Supplementary Fig. 1 and all data are presented in Supplementary Table 1.

Hydrocarbon emission
The predicted increase in monoterpene concentration is visible throughout the lower part of the troposphere and peaks during summer. The maximum increase occurs around noon at the surface layer, and later in the day in the upper levels of the boundary layer (Fig. 1b). During the spring and summer months, the strongest increase, with close to doubled concentrations, was Monoterpene increase by temperature. a The evolution of monoterpene concentrations in 2016 for the BASE case and the difference of the latter between the E6C6 scenario at the SMEAR II in Finland averaged for 3 km height above the surface. b The percentage change in vertical distribution of monoterpene concentration between scenario E6C6 and BASE relative to the maxima at each height level, averaged for daily distribution for the whole year 2016. The black line gives the yearly averaged planetary boundary layer height at SMEAR II. c The zonal mean values of monoterpene concentrations averaged vertically from 0-5 km (solid line) and 5-10 km (dashed line) for the Northern hemisphere calculated with the TM5 model for the BASE scenario. Added are the relative differences between the E6C6 and E0C6 scenario. The scenario abbreviations EXCX stands for E = emission, C = chemistry and X = temperature increase in the related module(s).
during daytimes, at altitudes above 2 km ( Supplementary Fig. 2b, c). However, the impact of monoterpenes on the tropospheric chemistry above the boundary layer during these seasons is small. This can be explained by the more than one order of magnitude lower monoterpene concentrations above 2 km compared to values in the mixed layer. This is due to the short atmospheric lifetime of monoterpenes with respect to reactions with OH, ozone and the nitrate radical (NO 3 ), which result in small fluxes out of the boundary layer ( Supplementary Fig. 1). In winter and autumn, the relative increase of the monoterpene concentration is temporally and spatially more even, with values around +40% to +50% for the E6C6 scenario when compared to BASE ( Supplementary Fig.  2a, d).
Next, when restricting the temperature increases to the emission module only (E2C0, E4C0 and E6C0 scenarios) and comparing them to the runs with increased temperatures in both the emission and chemistry modules (E2C2, E4C4 and E6C6 scenarios), the monoterpenes concentrations increased by 1.3-6% (Supplementary Table 1). This is due to two overlapping and contrary effects. First, the rates of the bimolecular reactions of OH with monoterpenes decrease with increasing temperature, which would lead to longer lifetimes of the terpenes. Second, monoterpene oxidation also leads to some production of OH (often called OH or HO x recycling or regeneration). This process is also temperature dependent, and in this case affects the terpene concentrations more. A more detailed discussion about these two mechanisms, and a case study for isoprene, is provided in the Supplementary Discussion.
Earlier research 24 found a global increase of 19% in monoterpene emission in the 2100 IPCC A1B climate scenario (1.8 K warming), while in another study 25 the global increase of monoterpene emission was 58% in the 2100 IPCC A2 climate scenario (4.8 K warming). In a recent publication 26 , Hantson et al. applied the dynamic global vegetation model LPJ-GUESS to calculate emission estimates for monoterpenes and isoprene over the period 1990-2100. Their results predict a 23% increase of global monoterpene emissions in the 21st century under the RCP 4.5 climate change scenario, with strong increases in monoterpene emission capacities towards the poles and decreases in the tropics. Based on their results, the E4C4 scenario can be the closest estimate for changes in emissions for the boreal area. However, in case that currently observed Arctic warming continues during the next decades, an increase of up to 6 K, corresponding to our E6C6 scenario, is definitely feasible in the boreal area and we applied the 6 K scenarios for the TM5 model simulations. Figure 1c presents the TM5 zonal mean values of monoterpene concentrations and the relative increase between the E6C6 and E0C6 scenarios averaged from ground to 5 km (solid lines) and 5-10 km (dashed lines) for the Northern hemisphere. The spatial and temporal distribution in the difference of monoterpene emissions for the boreal forest is provided under Supplementary  Fig. 3. The simulations predict an overall increase in Northern latitudes above 50°N of 29% and 42% in the first and second 5 km interval, respectively, and a significant rise of monoterpene concentrations in the high latitudes (>70°N) of 12% and 33% for the same height intervals. Although the absolute concentrations of the terpenes in this region are orders of magnitudes lower compared to the middle latitudes, small changes in the relatively clean Arctic atmosphere can have strong impacts on the chemistry in this area. Figure 2 presents the impact of the increased BVOCs concentrations on the hydroxyl radical distribution in the first 3.0 km of the troposphere calculated with SOSAA. Plots 2 A to 2D illustrate the percent change of the hydroxyl radical in the four scenarios where the temperature was increased in the emission and chemistry modules. The vertical distributions of the OH concentration for the BASE run averaged for daily concentrations for the four seasons in year 2016 are provided in Supplementary Fig. 4. The main loss of OH in all the scenarios occurs in the first 2 km above the ground during daytime. This is related to the relatively high monoterpene concentrations at this height interval. There is also a slight increase in the OH concentration during the night, which is related to the production of OH through the reaction of monoterpenes with ozone 27 . This is the main nighttime source for the hydroxyl radical. However, it only has a marginal impact on the diurnal OH budget, as the concentrations are about one order of magnitude smaller compared to daytime values. The overall OH differences relative to the BASE run averaged for the whole year 2016 at the selected height interval (0m-3km) are −1.2%, −2.4% and −3.5% for the scenarios E2C2, E4C4 and E6C6, respectively.

Impact of enhanced BVOC concentrations on OH concentrations
A closer look at the seasonal trends ( Supplementary Fig. 5) shows that the highest loss of OH in the boundary layer is observed during daytime in spring and summer. During this time, the change of the hydroxyl radical in the E6C6 simulation, relative to the BASE run, reaches values down to −30%, whereas in winter and autumn the model shows an overall increase of OH. However, the hydroxyl radical in the boreal region has the greatest impact on atmospheric chemistry during spring and summer months, as its main production mechanism involves O 3 photolysis, which requires solar irradiance. This leads to more than one order of magnitude higher OH concentrations during these periods compared to autumn and winter ( Supplementary Fig. 4). Hence, any impact through changed OH concentrations on the atmospheric chemistry will be most pronounced in spring and summer, which correlates strongly with the calculated increase of the monoterpene concentrations in the scenarios E2C2, E4C4 and E6C6.
By restricting the temperature increase to the emission module only (similar as discussed above for the monoterpene concentrations), or to the emission module and the first bimolecular reactions of OH, O 3 and NO 3 with the monoterpenes, the OH concentrations show an additional drop of about 1.3% when comparing E6C0 to the E6C6 scenario, and 0.7% when comparing E6CS to the E6C6 scenario (Supplementary Table 1). This shows that about half of the higher OH-loss in simulation E6C0 is related to ignoring the response of the OH-recycling to the increased temperatures for the monoterpenes currently included in MCM, while the other half comes from the higher rates of the initial bimolecular reactions. Thus, increasing the temperature also in the chemistry module buffers the decrease in the OH concentrations by approximately 37% (see detailed discussion in the Supplementary Discussion). An opposite OH trend is predicted when we only increase the temperature in the chemistry module, and let all other processes be equal to the BASE setup. The OH concentrations in the corresponding scenarios E0C2, E0C4 and E0C6 increase by 0.3, 0.6 and 0.8%, respectively.
In order to investigate the impact of the temperature increase on the hydroxyl radical concentrations in the Northern hemisphere (NH) we applied the TM5 model and calculated the relative drop of OH between the scenarios E6C6 and E0C6. Figure 3 shows the zonal averaged OH change for all seasons with a clear effect of the higher terpene concentrations at latitudes >50°N up to the tropopause. The mean decrease of OH for the whole year in the whole NH is −0.49% but doubles (−0.95%) if averaging only for latitudes north of 50°. One interesting feature is the influence of the terpene emissions on OH concentrations in the Arctic at higher altitudes during all seasons besides winter. As shown in Fig.  1c the terpene concentrations here are down to values below one million molecules per cm 3 , but still marginal rise of the absolute terpene concentrations decreases the OH by 1-2%. As this area is relative clean compared to the lower latitudes, a drop of OH by more than 1% can have a strong impact on the Arctic chemistry.
Further we explored the impact of the simplified chemistry scheme from TM5 on the OH budget by implementing the chemistry scheme from TM5 into SOSAA (named TM5-Chem). In the BASE simulations both models, SOSAA and TM5-Chem, show similar OH yearly mean concentrations between 0-3 km of 8.6 × 10 5 molecules cm −3 and 7.6 × 10 5 molecules cm −3 for TM5-Chem and SOSAA runs, respectively. Next, we compared the changes of OH concentration (relative to the BASE runs) when increasing the temperature only in the chemistry module, or in the chemistry and the emission modules of SOSAA and TM5-Chem ( Supplementary Fig. 6). In the two scenarios E0C6 and E6C6, TM5-Chem has much lower impact on OH concentrations than SOSAA. The SOSAA model with the chemistry scheme from MCM estimates a decrease of OH concentration by −3.5% when the temperature in the emission and the chemistry module is increased by 6 K (E6C6), whereas the TM5-Chem simulation only predicts a marginal OH drop of −0.11%. The difference between the two models is mainly related to the parameterized OHrecycling in the ozonolysis of monoterpenes in the TM5 chemistry scheme, which estimates higher OH concentrations during nighttime and a slight drop of OH during the day. When increasing the temperature only in the chemistry module (E0C6) the difference between SOSAA and TM5-Chem runs are smaller (OH change in TM5-Chem = 0.47% and in SOSAA = 0.84% for E0C6 runs). This strengthens our conclusion that the OH-recycling in the TM5 chemistry scheme is the main reason for the lower drop of OH concentrations in TM5-Chem. The OH concentration calculated by SOSAA is over seven times more sensitive to the temperature effects in the two scenarios than those calculated by TM5-Chem (E0C6-E6C6 for TM5-Chem = 0.58% and for SOSAA = 4.34%), highlighting a crucial difference between a simplified chemistry scheme applied in TM5 and a near-explicit chemistry scheme as MCM. This has a strong impact in future climate predictions for the OH concentrations when an increase in temperature is not considered in the estimation of the monoterpene emissions and could lead to an overestimation of the oxidation capacity of the lower atmosphere.

Methane concentration and lifetime increase
Although the percentage change of the OH concentration between the scenarios when the temperature increase is included in the emission module or not seems to be only a couple of percent compared to the high increase in monoterpene concentrations, it still has a substantial impact on the lifetime of methane. To estimate this impact, we applied the annual mean mixing ratio of CH 4 (1.95 ppm at SMEAR II 28 ), which is calculated from the average of three height levels (16.8 m, 67.2 m and 125 m above ground). We then applied the calculated hydroxyl radical concentrations from SOSAA for 2016 and calculated the time until the methane concentration reached 1/e of the original value, without considering the impact of increased temperature on other processes (e.g., meteorology or deposition). Here, the global atmospheric lifetime (also called turnover time) is the ratio of the number concentration of methane and the total rate of removal by OH from the reservoir 29 . The temperature used to calculate the CH 4 + OH reaction rate was then used in the associated scenario (i.e., the actual measured temperature at SMEAR II, raised by either 0, 2, 4 or 6 K). In this way, we are able to predict the lifetime of methane relative to the OH concentrations calculated in the different scenarios. The results are presented in Fig. 4, with an inline plot showing the seasonal percentage loss of methane averaged over all scenarios. Our calculated "local lifetime" is not intended as a prediction of the global atmospheric lifetime of CH 4 , but rather as an indicator of the local changes in the CH 4 oxidation rate.
The higher ambient temperatures in the scenarios EXCX (X = 2, 4 and 6) show a strong impact on the atmospheric lifetime of methane when compared to the scenarios where the temperatures were only increased in the chemistry module: E2C2 vs. E0C2~+3.4%, E4C4 vs. E0C4~+5.8% and E6C6 vs. E0C6~+11.4%. Another way of interpreting these changes is that the temperature increase tends to have a chemical influence of decreasing CH 4 lifetime, which is roughly offset by the emissions impact, such that across the simulations BASE, E2C2, E4C4, E6C6, the change in CH 4 lifetime is quite small.
To investigate the impact of these two scenarios (with and without considering temperature increase in the emission of terpenes) for the whole Northern Hemisphere we calculated the methane concentrations for E0C6 and E6C6 with the TM5 model. Figure 5 presents the zonal and yearly averaged relative methane concentrations change between the two simulations. The model predicts an overall increase of CH 4 concentration from ground to 10 km of 0.024% per year for the whole NH and 0.027% for the latitudes North of 50°. Taking into account that the yearly increase of atmospheric methane between 2008-2017 was about 0.35% (6.6 ± 0.3 ppb yr −1 ) 30 , there will be a significant underestimation of calculated methane growth rate by 7%, if the temperature effect is only considered in the chemistry but not in the emissions. This clearly shows that increased emissions of monoterpenes in the boreal region related to higher temperature in future climates will impact the methane concentration and lifetime in the whole Northern Hemisphere and not only in the high latitudes. A recent study on climate-driven chemistry in Earth System models 31 made similar conclusions, that the overall methane lifetime is expected to increase in a warmer climate due to increased BVOC emission.

Comparing SOSAA and TM5 results
The TM5 and SOSAA predicted changes of OH are in good agreement with similar patterns for the different seasons. However, when comparing the annual averaged differences in OH change at the location of the SMEAR II between the E6C0 and the E6C6 scenarios from SOSAA and TM5 we get a 56% lower value for the global chemistry model ( . The main reason for this difference is the impact of the two different applied chemistry schemes as all other modules and processes are unchanged. SOSAA is setup with nearly explicit chemistry from the Master Chemical Mechanism whereas TM5 is using a more simplified chemistry scheme for computational reason (see more details under methods).
Concerning how well the methane concentration or lifetime change prognosticated by TM5 and SOSAA between two scenarios are comparable requires a more detailed analysis (Fig.  6). SOSAA predicts an increase of methane lifetime of 11.4% when comparing the E6C6 and E0C6 scenarios, and TM5 a yearly increase of methane concentrations of 0.024% for the whole NH. If we consider the calculated lifetime of 8 years for methane in the E6C6 simulation by SOSAA and the 11.4% increase during this period, we get a 1.61% increase of methane concentrations per year at SMEAR II, the location for the SOSAA simulations (Provided in the Supplementary Methods is the theoretical proof for relating CH 4 concentration changes and to its lifetime). By further considering the percentage of the whole NH area that is boreal forest (4.7%) and only the outcome from TM5 in the first 3 km, we get a concentration change of methane for the Northern hemisphere in the first 3 km of 0.076 and 0.034% yr −1 by SOSAA and TM5, respectively. This reflects the same fraction of 44% that the two models predict in the OH drop. This further points to the conclusion that a more detailed chemistry scheme in TM5 than the one currently applied would lead to a more than double methane concentration change between the two scenarios.
To calculate the methane radiative forcing from our estimated yearly change in methane concentrations, we applied the formula from Etminan and co-workers 32 , which states that the change in radiative forcing of methane, ΔRF CH4 = 4.48 × 10 −4 Wm −2 ppb −1 . This leads to a ΔRF CH4 of 6.1 mWm −2 and 1.9 mWm −2 for the SOSAA and TM5 outcomes, respectively. Compared to the current methane radiative forcing 1 of 0.61 Wm −2 it points to an increase of 0.17% and 0.05% in ΔRF CH4 per year. Reflecting that future climate predictions are performed for decades or longer this could significantly impact future climate estimates if not considered.
Both models, although quite different in spatial and temporal resolution and in complexity of the applied chemistry and emission, lead to the similar result which is that future climate studies need to include interactive emissions and up-to-date chemistry schemes. This would ensure that the decrease of  methane concentrations in warmer climates through the temperature dependent chemical reactions will be compensated by higher terpene emissions and subsequent lower OH concentrations. In the case that a model doesn't take this into account, then the predicted radiative forcing will be underestimated for methane.

DISCUSSION
Compared to 1850-1900, global surface temperature averaged over 2081-2100 is very likely to be higher by 1.0°C to 1.8°C under the very low GHG emissions scenario considered, by 2.1°C to 3.5°C in the intermediate scenario (SSP2-4.5) and by 3.3°C to 5.7°C under the very high GHG emissions 33 . Moreover, it is nearly certain (99-100% probability) that the Arctic will continue to warm more than global surface temperature, with high confidence above two times the rate of global warming 33 . The feedbackincreased emission of BVOCs through climate warming in the boreal region, with a significant impact on the hydroxyl radical concentrations (decreased atmospheric oxidation capacity), followed by decreased methane oxidation-is only considered by selected highly coupled Earth System Models. Because the simulation of the feedback requires interactive biogenic VOC emissions, an interactive tropospheric chemistry scheme and an emission-driven methane cycle, it is likely that this complete feedback mechanism as such is missing from most Coupled Model Intercomparison Project Phase 6 (CMIP6) models 34 .
However, as pointed out in two recent publications 35,36 , a 4% reduction in OH, like predicted in our simulations comparing scenario E0C6 and E6C6 (4.3%), is equivalent to a 4% increase in methane emission, or 22 Tg/year (total annual global emissions of 548 Tg of CH 4 in the 2000s) 5 . We conclude that by ignoring the higher emission rates of monoterpenes from the boreal forest under future climate scenarios, will lead to underestimate the concentration, lifetime and climate impact of methane and thus to wrong conclusions about the future radiative forcing of methane. With potentially increasing anthropogenic emissions 37 together with predicted tree-line advance in the Arctic and subarctic regions 38,39 , and rise in methane emissions in the Arctic region through melting permafrost 40 , the effect of this feedback can become even stronger, and needs to be considered carefully in climate simulations.

METHODS
Model-SOSAA SOSAA (a model to Simulate the concentrations of Organic vapours, Sulphuric Acid and Aerosols) is a one-dimensional (1-D) chemical transport model which was first developed to investigate the BVOC emissions, the vertical transport and chemical reactions of trace gases within and above a boreal canopy in the planetary boundary layer (PBL) 20 . Since then, SOSAA has been applied in several cases to study, e.g., the oxidation of SO 2 by stabilized Criegee intermediate (sCI) radicals 41 , different monoterpene emission modules 42 , the aerosol dynamics in the PBL 43 , the characteristics of OH-reactivity 11,16,44 , the dry deposition processes of O 3 45 and BVOCs 46 . SOSAA has recently been validated against long-term measurements at SMEAR II for meteorological data and several selected atmospheric compounds, including monoterpenes and the hydroxy radical 47 . This manuscript also provides a sub-chapter discussing the uncertainties when applying SOSAA at the SMEAR II station.
SOSAA is written in Fortran and parallelized with MPI (Message Passing Interface). The source code is stored in an online git repository at https://bitbucket.org/ifbird/sosaa/src/master/, and continuously developed and maintained by the Multi-Scale Modelling Group at the University of Helsinki.
The current version of SOSAA contains five coupled modules which are described below: 1. The meteorology module is modified from a 1-D version of the three-dimensional (3-D) PBL meteorology model SCADIS (scalar distribution) 48 . This module is employed to calculate the prognostic equations for horizontal wind vector (u and v components), air temperature and absolute humidity. The upper boundary values of these prognostic variables are taken from the spatially and temporally interpolated ERA-Interim reanalysis data which are provided by the European Centre for Medium-Range Weather Forecast (ECMWF) 49 . In the lower part of the model domain from 4.2 m to 125 m above the ground, the wind vector, air temperature, and absolute humidity are nudged to the vertically interpolated measurement data at SMEAR II with a nudging factor of 0.05, which can represent the force of regional transport. At the canopy top, the long-wave radiation is provided by the ERA-Interim dataset. The incoming short-wave radiation is provided by the measurement data at SMEAR II. Then a radiative transfer module from the ADCHEM model 50 is used to split the observed short-wave radiation into a direct part and a diffuse part, both including a downward and an upward radiation component. This radiative transfer module uses the quadrature two-stream approximation scheme which was developed by Toon and co-workers 51 . The incoming downward direct and diffuse radiations, as well as the measured soil properties (soil temperature, soil water content and soil heat flux) are directly read in as input for the meteorology module. All of the meteorological input data mentioned above are linearly interpolated to 10 s time resolution to match the simulation time step.  45,46 . The Henry's law constants and reactivity factors of individual species in the chemistry scheme are calculated according to the methods described in a recent publication 46 . 5. The aerosol module is based on UHMA (University of Helsinki Multicomponent Aerosol model) 59 and implemented into SOSAA 43 . This module describes the aerosol dynamics, including the nucleation, condensation, coagulation, and deposition of aerosol particles. In this study, this module is not applied.
The model domain is usually set to 51 vertical logarithmically distributed layers, ranging from the surface up to 3 km to cover the whole PBL and a part of free troposphere. The simulation time step is 10 s for the meteorology module and 60 s for other modules. An implicit time integration method is applied to ensure the stability of computation. The output time step is 30 min.
Ozone concentrations applied in the model as mentioned under point 3 above were taken from the measurements at SMEAR II and not calculated. We note that using O 3 from input values probably underestimates the monoterpene increase, because in reality more monoterpenes emissions mean less O 3 (as the monoterpene + O 3 reactions consume it), meaning even higher monoterpene concentrations.

Model-EC-EARTH/TM5
TM5-MP (Tracer Model 5, Massively Parallel version) is an offline global chemical transport model, which is driven by the input meteorological and surface fields derived from the ECMWF (European Centre for Medium-range Weather Forecasts) ERA-Interim reanalysis datasets 60 . In this study, the modified version of chemistry scheme CB05 (carbon bond mechanism) 61 and the aerosol model M7 62 were used as in the default TM5-MP configuration.
The natural emissions of isoprene and monoterpene are prescribed as monthly mean values derived from MEGANv2.1 63,64 . A diurnal cycle is then applied to the monthly mean emissions to represent the effects of diurnal variations of radiation and temperature. The biomass burning emissions of monoterpenes and isoprene are applied from the datasets provided by van Marle 65 without diurnal cycle applied. The natural emissions of CH 4 are prescribed with the estimates for the year 2000 66 . The emissions of oceanic DMS (dimethylsulfide), mineral dust and sea salt, as well as NO x production by lightening are calculated online. The other natural emissions, such as terrestrial DMS emission, soil emissions of NO x and NH 3 , oceanic emissions of CO, emissions of NMVOCs (non-methane VOCs) and NH 3 , volcano emission of SO 2 , were prescribed as in van Noije 67 . The datasets provided by the Community Emissions Data System (CEDS) 68 and the BB4CMIP6 (biomass burning emissions for Coupled Model Intercomparison Project Phase 6) 65 are used for the prescribed anthropogenic and biomass burning emissions, respectively.
Specifically, for CH 4 , which is the focus of this study, its mixing ratio in the lower troposphere is nudged to the zonal mean surface mixing ratio provided by CMIP6 69,70 with the relaxation parameter of 2.9 day 71 . Here the lower troposphere refers to the domain from the surface to the highest full-level pressure layer above 550 hPa for a surface pressure of 984 hPa. In the stratosphere, CH 4 is nudged to an annual climatology value which is calculated as multiplying the measurement data from HALOE (Halogen Occultation Experiment) 72 by a scale factor. The scale factor is the ratio of the global mean surface mixing ratio between the annual mean value and the average value during 1991 to 2000. More details were described in van Noije et al. 73 .
In this study, we applied a horizontal resolution of 3 degrees in longitude and 2 degrees in latitude, and 34 hybrid-sigma levels in the vertical direction. The time step is one hour. All the simulations were run for the year 2009 with a one-year spin-up. TM5-MP was run in CSC (Finnish IT Center for Science) Puhti supercomputer. We applied 90 CPU cores for each run and the simulation time was about 10 h per simulation year.
In total, six simulation cases were conducted, which are BASE, E0C6, E6C6 and their non-nudging counterparts BASE-NN, E0C6-NN, E6C6-NN (Supplementary Table 2). Here non-nudging (NN) means the nudging of CH 4 to the zonal mean surface mixing ratios in the lower troposphere was switched off. The BASE case is the control case with the TM5-MP default configurations. In the four C6 cases, including E0C6, E6C6, E0C6-NN and E6C6-NN, the air temperature used in the chemistry mechanism was increased by 6 K over the whole globe, representing a warming future scenario. In the two E6 cases, including E6C6 and E6C6-NN, the monoterpene emission rate was increased by 71.6% throughout the boreal forest region as calculated from the SOSAA runs. The region was defined here as the area north of 50°N with the evergreen and deciduous needle-leaf trees covered during the whole year in 2009, which had similar vegetation cover as that at SMEAR II. Supplementary Fig. 5 shows the difference of monoterpene emissions over the boreal region between the E6C6 case and the BASE case in each month.

Uncertainties in EC-EARTH/TM5
The increasing rate of CH 4 mainly results from the loss of OH due to increased monoterpene emissions. First, considering the uncertainties of chemical reactions in the chemistry scheme, the annual mean tropospheric OH concentration uncertainty was estimated as 16% 74 . Secondly, although the monoterpene emissions vary in different large-scale models 75 , the emission rates calculated here are still the same in different models since it is a relative change. Therefore, the uncertainty range of the increasing rate is estimated in the range of 0.84 to 1.16 times of 0.024%, which is 0.020% to 0.028% and the corresponding uncertainty range of the increasing rate of radiative forcing lie in 0.011 to 0.015 Wm −2 per year. Besides, due to the limit of the simulation period, the uncertainty from interannual variability should also be investigated in future.

Station-SMEAR II
The measurements for 2016 applied in SOSAA were conducted at the SMEAR (Station for measuring Ecosystem-Atmosphere Relation) II station located in Hyytiälä (61°50'51''N, 24°17'41''E), Southern Finland 76 . The station is surrounded by 56 years old (in 2018) pine dominated forest, which also contains Norway Spruces and deciduous trees 54 . SMEAR II is a unique field station with continuous measurements of physical, chemical and biological phenomena, processes and interaction between these elements. Detailed description of the site can be found at the SMEAR II website (https://www.atm.helsinki.fi/SMEAR/index.php/ smear-ii).

Data processing
The global atmospheric lifetime is a general term used on various time scales for characterizing the rate of processes affecting the concentration of trace gases. In this study, we calculated the "Turnover time" (also called global atmospheric lifetime) 30 as the ratio of the average methane concentrations and the total rate of removal by the reaction with OH. The applied rate for the reaction of OH + CH 4 was k OH+CH4 = 1.85E-12 * EXP(-1690/T) (MCM version 3.3.1) 55 . The temperature applied here was from the 2016 BASE run, with a 2, 4 and 6 K increase as determined by the simulated scenario. In our study, we applied the OH concentrations obtained from the model for 2016 with a timestep of 30 min and calculated the time until the initial methane concentration reached a value <1/e of the original value. In this way, we obtained the local/ regional lifetime of methane for the different OH concentrations from our selected scenarios.

DATA AVAILABILITY
All data that support the findings of the study are available on request from Michael Boy (michael.boy@helsinki.fi).

CODE AVAILABILITY
The SOSAA code used for most of the model simulations in this study is available on request from Michael Boy (michael.boy@helsinki.fi). Access to TM5-MP can be obtained by emailing Philippe Le Sager (sager@knmi.nl).