Probabilistic projections of 21st century climate change over Northern Eurasia

We present probabilistic projections of 21st century climate change over Northern Eurasia using the Massachusetts Institute of Technology (MIT) Integrated Global System Model (IGSM), an integrated assessment model that couples an Earth system model of intermediate complexity with a two-dimensional zonal-mean atmosphere to a human activity model. Regional climate change is obtained by two downscaling methods: a dynamical downscaling, where the IGSM is linked to a three-dimensional atmospheric model, and a statistical downscaling, where a pattern scaling algorithm uses climate change patterns from 17 climate models. This framework allows for four major sources of uncertainty in future projections of regional climate change to be accounted for: emissions projections, climate system parameters (climate sensitivity, strength of aerosol forcing and ocean heat uptake rate), natural variability, and structural uncertainty. The results show that the choice of climate policy and the climate parameters are the largest drivers of uncertainty. We also find that different initial conditions lead to differences in patterns of change as large as when using different climate models. Finally, this analysis reveals the wide range of possible climate change over Northern Eurasia, emphasizing the need to consider these sources of uncertainty when modeling climate impacts over Northern Eurasia.


Introduction
Northern Eurasia accounts for 60% of the land area north of 40 • N and includes roughly 70% of the Earth's boreal forest and more than two-thirds of the Earth's permafrost [1]. As a result, the region is a major player in the global carbon budget. Over the past century, Northern Eurasia has experienced dramatic climate change, such as significant increases in temperature, growing season length, floods and Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. droughts [2,3]. These changes have large environmental and socioeconomic impacts including forest fires [4], permafrost thaw [5], extensive land-use change and water management projects [1]. Further climate change could lead to significant releases of greenhouse gas (carbon dioxide and methane) to the atmosphere caused by severe permafrost thaw, increasing forest fires, changes in lake and wetland dynamics and changes in land cover. This implies a potential positive feedback cycle. For this reason, it is imperative to quantify the full range of possible climate change over Northern Eurasia.
There is a large uncertainty in future projections of global climate change (see literature review in [6]) and regional climate change [7][8][9] associated with the uncertainty in, among others, internal variability [10][11][12], climate sensitivity [13][14][15], model uncertainty [16,17], and scenario uncertainty [18,19]. Other attempts to derive probabilistic forecasts of regional climate change include [20,21]. When it comes to climate change impacts over Northern Eurasia, recent studies include rising methane emissions [22], vegetation change [23,24], agroclimatic potential [25] and near-surface permafrost degradation [26]. However, these studies, along with many others focused on Northern Eurasia or other regions, generally rely on a small ensemble of climate simulations that does not cover the full range of uncertainty. In particular, such studies do not consider all the major sources of uncertainty in future projections of climate change, and thus are likely to underestimate the range of climate change and its impacts over the region.
In this study, we attempt to simulate possible future climate change over Northern Eurasia, by computing probabilistic projections of 21st century surface air temperature and precipitation changes and considering four major sources of uncertainty, namely: (i) uncertainty in the emissions projections, using different climate policies; (ii) uncertainty in the climate system response, represented by different values of climate parameters (climate sensitivity, strength of the aerosol forcing, ocean heat uptake rate); (iii) natural variability, obtained by initial condition perturbation; and (iv) structural uncertainty using different climate models. Our focus is the Northern Eurasian Earth Science Partnership Initiative (NEESPI) domain, which extends from 15 • E in the west to the Pacific coast in the east and from 40 • N in the south to the Arctic ocean coast in the north.

Modeling framework
This work uses the Massachusetts Institute of Technology (MIT) Integrated Global System Model (IGSM) [27,6], an integrated assessment model that couples an Earth System Model of Intermediate Complexity (EMIC), with a two-dimensional zonal-mean atmosphere, to a human activity model. The IGSM includes a representation of terrestrial water, energy, and ecosystem processes, global scale and urban chemistry including 33 chemical species, carbon and nitrogen cycle, thermodynamical sea ice, and ocean processes. The IGSM has been used in EMIC intercomparison exercises [28,29] as well as to perform probabilistic projections based on uncertainties in emissions and climate parameters [6,19]. In version 2.2, the IGSM uses a two-dimensional mixed layer anomaly diffusive ocean model. In version 2.3, the IGSM uses a three-dimensional dynamical ocean model based on the MIT ocean general circulation model [30,31]. Different versions of the ocean model exist with different values of the diapycnal diffusivity, which leads to different rates of ocean heat uptake. In the IGSM2.3, heat and freshwater fluxes are anomaly coupled in order to simulate a realistic ocean state. Observed wind stress from six-hourly National Centers for Environmental Prediction (NCEP) reanalysis [32] is used to more realistically capture surface wind forcing over the ocean. For any given model calendar year, a random calendar year of wind stress data is applied to the ocean in order to ensure that both short-term and interannual variability are represented in the ocean's surface forcing. Different random sampling can be applied to simulate different natural variability [33].
Regional climate change is then obtained from IGSM simulations using two downscaling methods. A dynamical downscaling method relies on the MIT IGSM-CAM framework [33] that links the IGSM version 2.3 to the National Center for Atmospheric Research (NCAR) Community Atmosphere Model (CAM) version 3.1 [34]. New modules were developed and implemented in CAM to allow climate parameters to be changed to match those of the IGSM. In particular, the climate sensitivity of CAM is changed using a cloud radiative adjustment method [35]. In the IGSM-CAM framework, CAM is driven by greenhouse gas concentrations and aerosol loading computed by the IGSM model, as well as by IGSM sea surface temperature (SST) anomalies. Because the IGSM-CAM relies on one single atmospheric model and because the cloud radiative adjustment method used to change the climate sensitivity does not provide enough difference in the patterns of regional climate change (unlike the perturbed physics approach), we explore the uncertainty in regional patterns of change using a pattern scaling approach based on projections from different climate models. This statistical downscaling method is based on a Taylor-expansion pattern scaling algorithm [17] that extends the latitudinal projections of the IGSM two-dimensional zonal-mean atmosphere by applying longitudinally resolved climate patterns from observations and from climate model projections from the Coupled Model Intercomparison Project phase 3 (CMIP3). The representation of feedbacks in the pattern scaling method is limited to the choice and the configuration of the CMIP3 models. This two-pronged approach simulates regional climate change at 2 • × 2.5 • resolution based on IGSM probabilistic projections. It has been used successfully in previous work on the United States [9].

Description of the simulations
In this study, we analyze two emissions scenarios corresponding to a median unconstrained emissions (UCE) scenario where no policy is implemented after 2012 and a stabilization scenario where greenhouse gases are stabilized at 550 ppm CO 2 (660 ppm CO 2 -equivalent) by 2100. The stabilization scenario corresponds to the level 2 stabilization (L2S) described in [36]. The UCE and L2S scenarios are similar to, respectively, the Representative Concentration Pathways RCP8.5 and RCP4.5 scenarios [18]. A more in-depth comparison with the RCP scenarios can be found in [33].
First, probability density functions of climate parameters (climate sensitivity, strength of the aerosol forcing, ocean heat uptake rate) are computed by running a large ensemble of the 20th century climate with the IGSM and comparing the output with observations, while accounting for errors in observation and natural variability [37]. Then, for each emissions scenario, a 400-member ensemble simulation with the IGSM2.2 is run with Latin hypercube sampling (LHS) of climate parameters from their probability density functions [6,19]. The uncertainty in the carbon cycle is also taken into account by varying the rate of carbon uptake by the ocean and terrestrial ecosystem. This approach is in line with the development of prior distributions used to run large ensemble of future climate projections, described in [38]. Pattern scaling is then applied to each IGSM2.2 ensemble member based on the patterns of climate change of 17 CMIP3 climate models, following [17]. The resulting meta-ensemble is viewed as a 'hybrid frequency distribution' (HFD) that integrates the uncertainty in the IGSM ensemble and in the regional patterns of climate change of different climate models. Each model is weighed equally, similarly to [16].
Additional simulations are conducted with the IGSM-CAM framework in order to complement the statistical downscaling with simulations using a three-dimensional atmospheric model. To limit the number of IGSM-CAM simulations, we use one particular version of the IGSM2.3 with an ocean heat uptake rate that lies between the mode and the median of the marginal posterior probability density function obtained with the IGSM in [13]. Then we choose three values of climate sensitivity (CS) that correspond to the 5th percentile (CS = 2.0 • C), median (CS = 2.5 • C), and 95th percentile (CS = 4.5 • C) of the marginal posterior probability density function with uniform prior (integrated over the net aerosol forcing). The values of climate sensitivity agree well with the conclusions of the Fourth Assessment Report (AR4) of the Intergovernmental Panel on Climate Change (IPCC), which finds that the climate sensitivity is likely to lie in the range of 2.0-4.5 • C [39]. The value of the net aerosol forcing is then chosen from the bivariate marginal posterior probability density function with uniform prior for the climate sensitivity-net aerosol forcing (CS-F aer ) parameter space [33], with the objective to provide the best agreement with the observed 20th century climate change. The values for the net aerosol forcing are −0.25 W/m 2 , −0.55 W/m 2 and −0.85 W/m 2 , respectively, for CS = 2.0 • C, CS = 2.5 • C, CS = 4.5 • . These three sets of climate parameters are shown [33] to reproduce the median, and the 5th and 95th percentiles of the probability distribution of 21st century global mean temperature change obtained in the IGSM ensembles previously mentioned [6]. We refer to these simulations are low, median and high IGSM-CAM simulations. Finally, 5-member ensembles were carried out for each choice of climate parameters and emissions scenarios using different initial conditions and random wind sampling (referred to as initial conditions in the remainder of the paper). Further details on the IGSM-CAM simulations used in this study can be found in [33].
In total, this study is based on 13 600 IGSM-HFD simulations and 30 IGSM-CAM simulations, providing an unprecedented ensemble of simulations using both dynamical and statistical downscaling. It should be noted that the probabilities of future climate projections presented in this study are dependent on the choice of climate model due to structural irreducible errors [38]. From here on, we refer to IGSM-HFD simulations corresponding to the 5th, median and 95th percentile of the NEESPI mean distribution of temperature or precipitation as, respectively, low, median and high IGSM-HFD simulations. Figure 1 shows 21st century time series of NEESPI mean surface air temperature and precipitation anomalies from present day for the IGSM-CAM and IGSM-HFD simulations. Even though the low, median and high simulations for each downscaling method are obtained from different distributions (NEESPI mean for IGSM-HFD and climate sensitivity for IGSM-CAM), the NEESPI means simulated by the two methods show a good agreement, especially for temperature. For precipitation, the IGSM-CAM tends to simulate stronger increases in precipitation than the IGSM-HFD simulations, most notably for the stabilization scenario. That is because the IGSM-HFD takes into account multiple models, some with lesser tendencies for increases in precipitation over Northern Eurasia than CAM. Overall, both downscaling methods show a large range of future warming (from 4.5 to 10.0 • C and from 2.0 to 4.0 • C for, respectively, the unconstrained and the stabilization scenario) and moistening (from 0.2 to 0.5 mm/day and from 0.05 to 0.25 mm/day for, respectively, the unconstrained and the stabilization scenario) over the NEESPI region. The stabilization scenario is always associated with a significant reduction in future climate change compared to the unconstrained emissions scenario. It should be noted that all of the IGSM-HFD simulations exhibit warming and moistening for both emissions scenarios, indicating the robustness of these tendencies amongst the CMIP3 climate models over the region. In addition, the IGSM-CAM simulations exhibit a much larger year-to-year variability than the IGSM-HFD, even in the mean of the 5-member ensemble based on different initial conditions. That is because the variability in the IGSM-HFD is solely driven by the IGSM two-dimensional atmosphere, thus underestimating local variability over the NEESPI region. The envelope of the 30 detrended IGSM-CAM simulations, which shows the unforced natural variability, shows a good agreement with the observed variability in NEESPI mean temperature and precipitation anomalies from 2000 to 2010. Finally, figure 1 reveals that the separation between the climate change and the unforced natural variability occurs at different times for temperature and precipitation, for each emissions scenario, and for the low, median and high simulations. In particular, the emergence of the anthropogenic signal from the noise occurs sooner for temperature than precipitation, and sooner for the unconstrained emissions than for the stabilization scenario. For the unconstrained emissions scenario, the warming emerges between 2020 and 2030 (based on both IGSM-CAM and IGSM-HFD low, median and high simulations) and the moistening between 2035 and 2055. For the stabilization scenario, the warming emerges between 2020 and 2040 and the moistening as early as 2030 and not quite yet by 2100 for the IGSM-HFD simulations. Another analysis comparing NEESPI mean changes in temperature and precipitation between the IGSM-CAM and IGSM-HFD is presented in figure 2. We compare IGSM-HFD frequency distributions of NEESPI mean temperature and precipitation changes for various periods of the 21st century with respect to present day to the range obtained from the IGSM-CAM simulations. Figure 2 further demonstrates the broad agreement between the two downscaling methods and the large range of plausible future warming and moistening over Northern Eurasia. A further analysis (not shown) reveals that the frequency distributions generally display a positive skewness and kurtosis (relative to the normal distribution). The positive skewness and kurtosis increase as the projections extend into the 21st century, and are larger for the unconstrained emissions scenario. The IGSM-CAM simulations also exhibit positive skewness, although it is more pronounced than for the IGSM-HFD. This can be explained by the fact that the IGSM-CAM simulations only consider one value of ocean heat uptake rate and that the marginal posterior probability density function with uniform prior for the climate sensitivity-net aerosol forcing (CS-F aer ) parameter space for this particular value of ocean heat uptake rate is itself skewed [33]. Figure 2 also illustrates the overestimation of precipitation increases from the IGSM-CAM compared to the IGSM-HFD. In addition, it shows that the full range of the IGSM-CAM simulations in the earlier part of the 21st century, largely driven by natural variability, can be as wide as the full range of the IGSM-HFD simulations. This suggests that the role of natural variability in driving the range of probable NEESPI regional changes is not negligible, especially for projections over the next few decades.

Results
The regional patterns of change over the NEESPI region simulated by the IGSM-CAM and IGSM-HFD approaches are then investigated. Figures 3 and 4 show maps of, respectively, 21st century changes in temperature and precipitation for the low, median and high simulations. Regional patterns of temperature changes agree well between the IGSM-CAM and IGSM-HFD, with the largest warming in the northern parts of the NEESPI region. For precipitation, there is also a broad agreement in the pattern of drying/moistening between the two downscaling approaches, with some drying in Eastern Europe and the southern parts of the NEESPI region and moistening in the northern parts. The IGSM-CAM simulations show similar patterns of temperature and precipitation changes, with larger magnitudes for higher climate sensitivities and emissions. This is because the IGSM-CAM relies on a single atmospheric model and because figures 3 and 4 show the average over the five initial conditions. Averaging over the different initial conditions filters out most of the natural variability, leaving only the human induced climate response, which displays similar patterns of change even with different values of climate sensitivity [35]. On the other hand, the IGSM-HFD simulations show larger differences in the patterns of change because they consider multiple models and thus include structural uncertainty. Figure 5 shows the impact of the initial conditions within the IGSM-CAM framework. Maps of 21st century changes in temperature and precipitation for the median simulation under the stabilization scenario and for different initial conditions reveal the significant role of natural variability in future climate projections over Northern Eurasia. With different initial conditions, the simulations show similar magnitudes in temperature and precipitation changes but very different patterns. The location of the maximum warming can differ significantly, from European Russia (initial condition 3) to Eastern Siberia (initial condition 5). Precipitation patterns are also strongly influenced by the initial conditions, with a significantly different extent of the drying pattern found over Eastern Europe and the southern parts of the NEESPI region. The location of the maximum moistening can vary widely, from Scandinavia (initial condition 4) to Northern China (initial condition 2). The impact of the model pattern in the IGSM-HFD approach is analyzed by plotting the median simulation under the stabilization scenario and the four surrounding simulations, corresponding to the 50.02th, 50.01th, 49.99th and 49.98th percentiles of the NEESPI mean distribution ( figure 6). The NEESPI mean of these five simulations is virtually identical and each simulation could be considered as the median. However, the associated pattern of change is often very different because the corresponding model used in the pattern scaling method is different. This leads to differences in temperature patterns similar to the initial condition analysis, with different locations of the maximum warming. For precipitation changes, the five IGSM-HFD simulations show less discrepancies than the    initial condition analysis, largely because three out of the five simulations rely on the same model, and because the other two are based on models that seem to have similar patterns of precipitation changes over Northern Eurasia. This is a surprising result that shows that the uncertainty in regional climate change simulated by ensembles based on initial condition perturbation and multimodel ensembles seems to compare well over Northern Eurasia.

Summary and conclusion
In this study, we present probabilistic projections of climate change over Northern Eurasia (NEESPI region) using the MIT IGSM downscaled via both a dynamical method (the IGSM-CAM framework) and a statistical method (pattern scaling). The analysis of the very large ensemble of simulations (a total of 13 630 simulations) shows that the uncertainty in the choice of policy and in the climate response (climate sensitivity, strength of the aerosol forcing and ocean heat uptake rate) results in a wide range of probable outcomes. It further shows that simulations with different initial conditions can lead to different patterns of change (even in the 20-year mean changes), as different as using different models. This is especially true for lower values of climate sensitivity and emissions scenarios with stringent stabilization of greenhouse gases. In addition, the precipitation change signal for the low simulation and stabilization scenario has not emerged from the noise even by 2100. The role of the uncertainty in natural variability shown in this study is in agreement with [11] that shows that natural variability contributes substantially to the uncertainty in climate change projections. This result also suggests that, for simulations with a relatively small warming (low climate response and small greenhouse gas radiative forcing), an ensemble based on initial condition perturbation could potentially be used within a single model as a substitute for a multimodel ensemble, even for end-of-century projections. However, this study suggests that at the scale of Northern Eurasia, the choice of policy is the largest source of uncertainty, followed by the climate parameters. This is in agreement with the findings of [9] for the United States. This findings is especially true for long-term projections that extend past 2050. Generally, the temporal changes in the contributions of the various sources of uncertainty are consistent with the works of [7] and [8].
It should be mentioned that this study suffers from limitations, such as the relatively low resolution of the IGSM-CAM and IGSM-HFD simulations or the absence of possible feedbacks (e.g., land-use change, aerosol-cloud interactions. . . ). Also, not all sources of uncertainty are considered, such as the uncertainty arising from different model resolution or the uncertainty involved in the pattern scaling method itself. Nevertheless, in light of these projections, it appears obvious that Northern Eurasia is at risk of substantial climate warming if mitigation policies are not implemented. Based on recent observed trends, such warming could lead to further widespread permafrost degradation and more intense and frequent forest fires [4], and potentially result in the release of large amounts of carbon and methane [40]. The simulations with different emissions scenarios, values of climate parameters, initial conditions and models show consistent patterns of drying in the southern parts of the NEESPI region, especially over Eastern Europe, and moistening over the rest of the region. These pronounced features indicate potential predictability in future precipitation changes over the region.
Overall, we recommend that when investigating climate change impacts over Northern Eurasia, studies consider at least the four sources of uncertainty analyzed in this study, namely: (i) uncertainty in the emissions projections, using different climate policies; (ii) uncertainty in the climate system response, represented by different values of climate parameters; (iii) natural variability, using different initial conditions; and (iv) structural uncertainty using different climate models. Furthermore, we suggest that probabilistic projections be used to drive impact models, even though we realize it would require large computing capabilities and would put a larger burden on impact modeling groups. Nonetheless, in light of this study, it appears evident that uncertainty in regional climate change projections is still large and should be accounted for systematically when estimating regional climate impacts. Because uncertainty in future climate projections is conditional on the methodological choice to derive probabilistic distributions and is affected by the model used, any impact analysis should explore the sensitivity of its results to different methodologies and models.