Large-scale emulation of spatio-temporal variation in temperature under climate change

Future temperature variations under greenhouse gas (GHG) emission scenarios are critical to assess possible impacts on human society and make reasonable mitigation policies. Due to the huge running cost, Earth system models (ESMs) may be difficult to flexibly provide the temperature projections following some specific emission pathways for empirical analysis. This study develops the mean and variability filed emulators in the high-resolution land grids to approximate the temperature behavior conditioned on GHG emissions in ESM. The emulator of mean temperature response is modeled as a function of GHG emissions to represent the expected values for ESM output, and the associated high-dimensional spatial dependence across grid points is estimated by the nearest-neighbor Gaussian process. The variability emulator is constructed with the residuals between the mean temperature response and the ESM output, and the associated space-time correlation structure is decomposed by principal component analysis and discrete Fourier transform. The analysis shows that the emulators trained with the runs of ESM only from part of representative concentration pathways can efficiently reproduce the temperature variations under different emission scenarios. The emulated gridded temperatures would be easily taken for climate impact and risk assessment, and be incorporated in the integrated assessment model for climate policy analysis.


Introduction
Future temperature variations under greenhouse gas (GHG) emission scenarios are concerned by climate change community to assess possible impacts on human society and make reasonable mitigation policies (Labriet et al 2015, James et al 2017, Agliardi et al 2019. At present, Earth system model (ESM) is state-of-the-art tool for revealing the linkage between temperature and GHG emissions. However, its huge running demands limit the use in the analysis of climate change impact and policy, and thus there is an increasing need for computationally cheap emulator to approximate ESM output (Link et al 2019, Beusch et al 2020. The emulated temperature response to GHG emissions, representing the current best understanding of the complex dynamical feedbacks of climate system, essentially provides scientific basis for estimating the likely impacts of climate change and choosing the proper emission reduction plan in the integrated assessment model (IAM) (Holden et al 2014, Hartin et al 2015, Foley et al 2016, Tebaldi and Knutti 2018, Tebaldi et al 2020. Given a specific emission scenario, the emulated range of temperature variations further facilitates quantifying the risks of climate change. Recent studies suggest that the statistical emulators are useful for efficiently reproducing the temperature behavior under emission scenarios (Castruccio et al 2019, Holden et al 2019). For example, Miftakhova et al (2020) employ an econometric model to map emissions to global mean temperature which is a crucial indicator of climate change, and such a low dimensional relationship can be directly implemented for investigating the dynamic interactions between climate and economy in some IAMs. In reality, however, the local climate information is more important for calculating the impacts of particular emission scenarios. The temperature transformation from global level to local level is commonly made by pattern scaling which assumes a linear relationship between global mean temperature and local temperature (Tebaldi and Arblaster 2014, Alexeeff et al 2018, Zelazowski et al 2018. Yet, this assumption is inconsistent with the fact that the temperature responses to GHG emissions are spatially varying over time (Holden and Edwards 2010). Therefore, an alternative is to model regional temperature directly conditioned on GHG emissions rather than global mean temperature. For instance, Castruccio et al (2014) employ an infinite distributed lag model to capture the nonlinear evolution of regional climate under the past trajectory of atmospheric carbon dioxide (CO 2 ) concentrations. Since temperatures correlate across locations, it is supposed to improve the statistical emulator's performance by including spatial dependence (Bao et al 2016).
There has been a considerable literature in spatial modeling for climate (Castruccio 2016, Lu et al 2018, Arisido et al 2019, Mukhopadhyay et al 2019, Bakar 2020. Usually, the dependence between local temperatures is captured by spatial covariance in a hierarchical framework (Gelfand et al 2005). Fitting such model, however, can be prohibitive with an increase in locations (i.e. the big N problem) (Eidsvik et al 2012). It means that the emulation of gridded temperatures considering spatial dependence is faced with computational intractability. A variety of solutions to this problem have been proposed. For example, the low rank methods, such as fixed rank kriging (Cressie and Johannesson 2008) and predictive processes (Banerjee et al 2008, Sahu and Bakar 2012), create approximation in a subset of locations to reduce the dimension of covariance matrix. In order to ensure adequate approximation, however, the computational cost becomes higher as the number of selected locations increases. Besides, the strong correlations between temperatures in neighboring locations could make these models perform poorly (Stein 2014). Alternatively, the sparse methods introduce sparsity in the covariance and precision matrices for efficient computation (Furrer et al 2006, Lindgren et al 2011, Liang et al 2013, Datta et al 2016, Guinness and Fuentes 2017, Katzfuss 2017. Heaton et al (2019) compare different approaches to show that the nearest-neighbor Gaussian process (NNGP) models yield highly competitive predictive performance and computation time. This method, providing fully model-based inference, could help emulate the temperature variations across grid points which consist of large-scale space-time processes.
This study intends to construct the emulators of temperature response to GHG emissions in the highresolution land grids, including both the mean and variability fields. First, we develop the emulator of mean temperature response as a function of GHG emissions using the NNGP to represent the expected values for ESM output. The departures of real ESM output from the mean response are viewed as variability. Thus, the residual temperature variability is then generated by the space-time decomposition based emulation technology (Link et al 2019). The simulated temperatures of ESM under representative concentration pathways (RCPs) are used for calibrating and verifying the emulators. The analysis shows that the emulator trained with the runs only from part of RCPs can efficiently reproduce the temperature variations under different GHG emission scenarios. The emulated gridded temperatures would be easily taken for climate impact and risk assessment, and be incorporated in the IAM for climate policy analysis.
The rest of this paper is organized as follows. Section 2 describes data sets and the models for spatio-temporal temperature emulation in our analysis. In section 3, we provide the emulated mean temperature response and residual temperature variability at the grid level. On the basis, the synthetic temperature fields are presented and the performance of emulators is evaluated. Finally, section 4 concludes the paper.

Data description
This study intends to emulate the annual temperature behavior with changing GHG emissions using the runs of ESM. As an emulation example, the single run under each RCP from Beijing Climate Center Climate System Model version 1.1 (BCC-CSM1.1) (Xin et al 2013) is used to conduct our analysis. Specifically, the simulated temperatures over the historical period  and under RCP2.6 andRCP8.5 (2006-2100) are taken to train the emulators. On the basis, the temperatures under RCP4.5 and RCP6.0 (2006-2100) are employed to evaluate the performance of emulators. Considering the demands for highresolution climate information, the annual temperatures from BCC-CSM1.1 are bilinearly interpolated onto the 0.5 • × 0.5 • land grids for emulation. The grid points cover eight regions (figure 1) which are modified according to the definition by Ruosteenoja et al (2003). In this study, the global mean refers to the average across all land grid points. As the driving factor of temperature variations, the annual global GHG emission (carbon dioxide equivalence) concentrations are taken from the RCP Database. 4

Emulator construction
In the next, matrices and vectors are written in bold letters, and vectors are assumed to be column vectors. A superscript T represents the transpose of a matrix or a vector.

Mean temperature response
The mean temperature is emulated as a function of carbon dioxide equivalence concentrations using a hierarchical Bayesian approach. We define CC(t) as the ratio of current (the tth year) and preindustrial (the year of 1765) carbon dioxide equivalence concentrations. For the given location of s i and year t, the short-and long-term carbon dioxide equivalence concentration covariates ( where w(k) = p k−2 (1 − p). The long-term carbon dioxide equivalence concentration covariate is assumed to combine a set of carbon dioxide equivalence concentrations in the prior years with an exponential decay weight. The parameter p captures the decreasing impacts of past carbon dioxide equivalence concentrations on current temperatures. As a result, the mean temperature response to the shortand long-term carbon dioxide equivalence concentrations in location s i and year t is expressed as where α(s i ) is spatial random effect, indicating local adjustment to the mean temperature. β 1 (s i ) and β 2 (s i ) are the regression coefficients associated with the short-and long-term effects of carbon dioxide equivalence concentrations. ε(s i ,t) is a sequence of identical and independently Gaussian distributed random errors with mean 0 and variance τ 2 . The spatial effects α(s)=(α(s 1 ), α(s 2 ), . . . , α(s n )) T with location vector s = (s 1 , s 2 , . . . , s n ) T are commonly specified as a Gaussian process with mean u and covariance matrix C, denoted α ∼ GP(u, C(·, ·; σ, φ)). The (i, j)th entry of exponential covariance matrix is Here, φ is decay parameter and ∥ s i − s j ∥is the great circle distance between the locations s i and s j . The effective range of ∥ s i − s j ∥ can be found in the work of Banerjee et al (2008). Since the large number of locations makes it difficult to estimate the covariance matrix, the NNGP (Vecchia 1988, Stein et al 2004, Datta et al 2016 is employed to make an approximation of C by generating the sparse covariance matrix. To assess the spread of covariate effects across locations, the regression coefficients with spatial structure are also specified as Gaussian processes, denoted β 1 ∼ GP(u 1 , C(·, ·; σ 1 , φ 1 )) and β 2 ∼ GP(u 2 , C(·, ·; σ 2 , φ 2 )). The spatially varying coefficients imply that the effects of carbon dioxide equivalence concentrations are supposed to vary smoothly in space. In addition, the regression coefficients without spatial structure are further specified as normal distributions, denoted β 1 ∼ N(µ 1 , Σ 1 ) and β 2 ∼ N(µ 2 , Σ 2 ). If the estimated variances of Σ 1 and Σ 2 are large, then effectively it indicates that each location is regressed independently. By contrast, small variances imply homogeneous responses to carbon dioxide equivalence concentrations. We use noninformative priors for all the parameters and employ Markov Chain Monte Carlo (MCMC) sampling to estimate posterior distributions. The convergence of the MCMC chain is evaluated by the potential scale reduction factor.

Residual temperature variability
The residual temperature is obtained by subtracting mean response from ESM output, implying the variation around the mean response. A crucial task for emulating such temperature variability is to generate the same distribution and space-time correlation structure as the residual temperatures. Following Link et al (2019), we use principal component analysis to decompose the residual temperatures r(t) = (r(s 1 , t), r(s 2 , t), . . . , r(s n , t)) T into wherer Ther(q) is empirical orthogonal function (EOF) computed by singular value decomposition, and it reveals the space structure uncorrelated over time. The λ(q, t) is the projection coefficient for EOF, and it can be generated independently according to equation (7). Specifically, we compute the discrete Fourier transform (DFT) of projection coefficients λ(q) = (λ(q, 1), λ(q, 2), . . . , λ(q, t)), i.e. Λ(f) = F(λ(q)). Then, the new Fourier transform Λ * (f) such that |Λ * (f)| = |Λ(f)| is generated by choosing its phases randomly on the interval [0, 2π]. According to the Wiener-Khinchin theorem, the reconstructed λ * (q) that is the inverse DFT of Λ * (f) has the same autocorrelation as λ(q). As a result, the residual temperature variability can be reproduced by equation (5) with λ * (q).

Mean field emulator
We separately establish the mean response emulators for the eight regions defined in this study regarding to the different climate characteristics. The emulators are trained with the RCPs providing possible high and low GHG emissions, so that it could be helpful for capturing the temperature variations within the range of these emissions. As suggested by Bao et al (2016), the parameter p for exponential decay weight is assumed to be greater than 0.95 to avoid the multicollinearity between covariates. In addition, the maximum effective range is considered as much as 50 • between two grid points, and the number of neighbors m = 10 is selected for estimation. We also increase the neighbors up to 20, but the associated verification performance is not improved. The medians of posterior distributions of model parameters are presented and used for emulation.
There are significant differences in the estimated parameters for the mean temperature response emulators across the regions. Table 1 shows the parameters for the emulators with spatially varying effects of carbon dioxide equivalence concentrations. It can be found that most of the estimated decay parameters are close to the assumed lower bound for the regions, which implies that temperatures correlate over a wide area. In ARL, NAM, NAS and SEA, the larger mean values u 1 suggest that the temperatures are sensitive to the short-term carbon dioxide equivalence concentrations. In addition, the larger variances σ 1 and σ 2 indicate the distinctly varying mean temperature responses to carbon dioxide equivalence concentrations. By comparison, the smaller spatial variability in the effects of carbon dioxide equivalence concentrations on the mean temperature response is found in AFR and PAC whose temperatures are generally higher and tend to have stronger correlations. Besides, the σ suggests that the baseline temperatures in NAM and SEA, which cover North America and South and East Asia, have higher spatial variability.
The assumption on the spatial structure in regression coefficients β 1 and β 2 is also relaxed. To account for the common information shared by the grid points within a region, the regression coefficients are modeled with normal distributions. As shown In addition, ARL also highly responds to the long-term carbon dioxide equivalence concentrations. Although the spatial dependence of regression coefficients is not considered, the heterogeneous temperature responses to carbon dioxide equivalence concentrations are suggested by the estimated variances of regression coefficients. Yet, the smaller variances Σ 2 in NAS and SEA imply the regional homogeneous responses to the long-term GHG emissions. For the emulators with and without spatially varying regression coefficients, the spatial effects α have similar patterns. The trained emulators for mean temperature response are used to generate the gridded temperatures under RCP4.5 and RCP6.0 over 2006-2100, and the associated emulation performance for each region is evaluated by the root mean square error (RMSE) and the mean absolute error (MAE) (table 3). Since the temperatures highly fluctuate over time, ARL, NAM and NAS have larger RMSE and MAE. By comparison, the temperature fluctuations are of narrow amplitudes in AFR, SAM and PAC, and thus the fitting errors are smaller. Generally, the emulators with and without spatially varying regression coefficients display similar fitting performance, and in many cases dropping the spatial structure in regression coefficients seems better. It implies that the varying coefficients can capture the heterogeneous temperature response to carbon dioxide equivalence concentrations at the grid level, and modeling their spatial dependence may not always be necessary. Therefore, the emulators with spatially varying intercepts only are used in our analysis. We compare the grid-level mean field temperatures with ESM temperatures during different periods (figure 2). It shows that the emulators perform temperature changes well over 2061-2100, and the correlations between the emulated and real temperature changes during 2081-2100 are 0.96 and 0.97 under RCP4.5 and RCP6.0, respectively. For the global level, the mean field emulators successfully capture temperature trends, and there are small differences in the average changes in the mean field temperatures and ESM temperatures (figure 3). The emulated global temperatures are lower than ESM temperatures before 2060 but higher after 2060 under two RCPs. Generally, the mean temperature emulators can be used to indicate global temperature variations that are crucial for the economic analysis in some IAMs.
The average changes in the mean field temperatures over 2006-2100 are displayed in figure 4. The overall spatial patterns of ESM temperatures are well captured by the mean field emulators. We use the absolute differences between the mean field temperatures and ESM temperatures to quantify the emulators' performance. The absolute temperature difference is less than 0.1 • C for 66.1% of the grid points  and less than 0.5 • C for 99.0% of the grid points, and the average for all the grid points is 0.1 • C. The larger differences are mainly in northeastern Canada and southern Greenland under RCP4.5, and in central Asia and eastern Russia under RCP6.0. On the whole, the mean temperature emulators can reproduce the spatial variations of gridded temperatures.

Variability field emulator
According to the emulated mean temperatures, the associated residuals are taken as temperature variability. An emulator of variability field is constructed with all gridded residuals. The spatial patterns of leading EOF modes of variability field are shown in figure 5. The four EOF modes explain 20.9%, 12.0%, 7.5% and  6.6% of total variance, respectively, and are characterized by the positive anomalies mainly in the northern hemisphere. For example, the anomalous positive values appear in southern Greenland and southern the United States in the first EOF (EOF1), and then expand to most of the United States, Europe and northern Asia in the second EOF (EOF2). However, there seems to be opposite patterns in North America  and most of Europe between EOF2 and EOF3. Also, the patterns of positive anomalies in Europe and most of Asia are generally opposite between EOF3 and EOF4. By keeping full set of EOFs and randomly modeling the projection coefficients of EOFs, the temperature variability field is emulated.
The performance in emulating temperature variability is evaluated by checking the associated statistical properties. We generate an ensemble of emulated residual temperature fields, and compare the standard deviations of emulated and ESM residual temperature fields at grid points (figure 6). It shows that the emulator can generally replicate temperature variability, and the ensemble of 50 emulations has better performance. As indicated by Link et al (2019), the statistical test suggests that there are very slight discrepancies between the emulated and real temperature variability fields, and the spatial and time correlations are retained in the emulations. Therefore, such emulator can reliably reproduce temperature variability, and we use it to emulate variability fields conditioned on the mean temperature fields.

Temperature emulation
The synthetic temperature variations at grid points are generated by the mean and variability field emulators, and the global mean temperatures under RCPs are presented in figure 7. It shows that the averaged emulations accurately reflect the temperature trends, and the emulation ranges well cover the ESM temperature variations. There are 95.8% and 99.0% of the temperatures within the emulation ranges under RCP4.5 and RCP6.0, respectively. Previous physical studies suggest that temperature change is highly related to the ratio between current and preindustrial GHG emission concentrations (Forster et al 2007, Millar et al 2017. Therefore, the present and historical changes in GHG emission concentration are taken as the drivers of mean field temperature change in this study. The emulators are trained with the ESM output in the historical period and under RCP2.6 and RCP8.5, so that it can reflect a wide range of mean field temperature changes with the associated GHG emission concentrations. This enables the mean field emulator to correctly perform the temperature trends under RCP4.5 and RCP6.0. The trained variability field emulator is the approximate reflection of temperature variations under RCP4.5 and RCP6.0, but shows reasonable coverage of temperature variations. As a result, the temperature emulators trained with the information from RCP2.6 and RCP8.5 basically reveal the temperature behavior under RCP4.5 and RCP6.0. It is also reliable to produce temperature variations under other emission scenarios using the trained emulators. Due to the cheap computation, it is easy to obtain a great amount of temperature runs. The generated temperature emulations would facilitate quantifying the climate uncertainty which is of great importance in the impact analysis.
For each emulation, the mean value and standard deviation of temperatures during 2081-2100 are calculated. Then, we average the mean values and standard deviations of 50 emulations to compare with those of ESM, and the corresponding differences are shown in figures 8 and 9. There are slight differences for most of the areas. The absolute difference in mean values is less than 0.5 • C for 96.0%/93.8% of the grid points under RCP4.5/RCP6.0, and that in standard deviations is less than 0.2 • C for 86.7%/92.8% of the grid points under RCP4.5/RCP6.0. The significant differences in mean temperature mainly appear in the United States and eastern Canada under RCP4.5, and northern Europe and Russia under RCP6.0. Overall, the emulated gridded temperatures well represent the ESM output.

Conclusions
The temperature response to GHG emissions is critical for the analysis of climate change impact and policy. This study constructs the mean and variability filed emulators in the high-resolution land grids to approximate the temperature behavior conditioned on GHG emissions in ESM. The mean temperature response represents the expected values for ESM output, and the associated spatial dependence across grid points is estimated by the nearest-neighbor Gaussian process model. The departures of real ESM output from the mean response indicate temperature variability, and the associated space-time correlation structure is decomposed by principal component analysis and discrete Fourier transform. Overall, the largescale temperature variations under climate change can be efficiently and accurately reproduced.
The emulators would bridge ESM and IAM, which is meaningful for analyzing the dynamic feedbacks between climate and economy. We train the emulators with ESM output to represent the current understanding of climate system. Therefore, incorporating the designed relationship between temperature and GHG emissions would bring the knowledge of complex processes from ESM into IAM. This provides solid temperature response to GHG emissions and helps quantify the high-resolution climate impacts, which makes the integrated assessment reliable. Besides, the emulation range characterizes the climate uncertainty. The decision makers would be richly informed by the range of temperature changes so that the robust policies could be obtained.
The extension of emulation could be easily made. The seasonal and monthly temperatures would be generated by including the seasonal and monthly components in the mean field emulator (Castruccio et al 2019). For other climate variable such as precipitation, the associated mean response function is supposed to include some additional driving factors (Castruccio et al 2014). In addition, the emulators are also appliable for other ESMs, and thus an ensemble of emulated temperatures would inform uncertainties from different ESMs. As suggested by Link et al (2019), a time varying scale factor may be needed to characterize the changing statistical properties in the variability field.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.