Simulation of future climate under changing temporal covariance structures

A growing body of evidence indicates that anthropogenic greenhouse gases are changing Earth’s climate, and that those changes may involve not only changes in climatic means but also in variability. Climate models may be informative about these future changes, but their use is complicated by the fact that they do not capture variability in current climate well. Many methods have therefore been developed to combine models and data in simulations of future climate, but current methods generally account only for changes in marginal variation and do not capture projected changes in correlation (spatial, temporal, spatiotemporal). We develop here a procedure to simulate future daily mean temperature that modifies climate observations based on changes in the mean and spectral density suggested by climate model output, and illustrate our methodology with projections from the CCSM3 (Community Climate System 3) climate model. We are able to simulate a future climate with changing temporal covariance while largely retaining non-Gaussian features of the observations. Our results suggest that in CCSM3, at most locations and most timescales, variability in daily mean temperature decreases under anthropogenic warming. The methodology presented here applies only to fully equilibrated future climate states, but may be extended to simulating transient states as well.


Introduction
With mounting evidence indicating that Earth's climate is changing (IPCC, 2007, and references therein), it is becoming increasingly important to understand the potential impacts of climate change on society.Impacts assessment requires projections of future climate under increased concentrations of greenhouse gases (GHGs).For example, understanding climate effects on food supply would require simulations of future temperature and precipitation for use in agricultural yield models.Crop yields, however, are highly nonlinear with temperature and precipitation and therefore are sensitive not only to climatological means but also to shortterm extremes (e.g., Schlenker and Roberts, 2009;Wheeler et al., 2000).In this context, climate must be understood as an underlying, multivariate, spatiotemporal probability distribution, for which weather is a random realization.Human societies can be impacted by changes of not only the mean, but of the entire probability distribution.
Changes in variability in both temperature and precipitation are physically plausible.For precipitation, standard physics would suggest increases in both spatial and temporal variability, with dry areas drier, wet areas wetter, and rainfall occurring in more intense events (Held and Soden, 2006).Changes in temperature variability are less robustly predicted.Some empirical studies suggest that temperature variability may already be changing in particular contexts, though some studies argue for increases and others for decreases (Karl et al., 1995;Timmermann et al., 1999;Schär et al., 2004;Hansen et al., 2012;Huntingford et al., 2013).For instance, in large portions of North America, a weakened polar jet stream (in part driven by reduced latitudinal temperature gradients) may (Francis and Vavrus, 2012) or may not (Barnes, 2013;Screen and Simmonds, 2013) lead to prolonged weather patterns that shift variation to lower frequencies.Studies of variability changes in climate models remain limited, and it is not yet clear which if any predictions are robust across models (see, e.g., Barnes and Polvani, 2013) or physical parameterizations (see, e.g., Hawkins et al., 2013).
The area remains one of active research.
One complication to analyses of potential future changes in climate variability is that while the deterministic climate models used for long-term climate forecasts appear to capture trends, they do not accurately reproduce observed current climate.These models, known as atmosphere-ocean general circulation models (AOGCMs), are physically based numerical simulations of transport of energy and moisture in the atmosphere and ocean, typically with separate submodels for the atmosphere, ocean, sea ice, and vegetation.Many AOGCMs successfully reproduce observed large-scale circulation, atmospheric structure, latitudinal temperature gradients, storm tracks, and quasi-periodic interdecadal phenomena such as the El Niño-Southern Oscillation.When driven with historical records of CO 2 and aerosol emissions due to human and volcanic activity, they also reproduce well the observed temperature trend of the last 2 centuries.Figure 1 demonstrates this ability to capture trends in the widely used Community Climate System 3 (CCSM3) model (Collins et al., 2006), which we use in examples throughout this manuscript.(See Appendix A for description of model and experimental runs, as well as observational data used in comparisons.)CCSM3 and other AOGCMs do not, however, perfectly reproduce either the mean or distribution of current climate.Model present-day global mean temperature (GMT) can be offset by several degrees from observations (again, see Fig. 1) and probability distributions of temperature and precipitation at individual locations do not match those of weather observations (Fig. 3, which shows marginal distributions in CCSM3 temperature output and observations for three representative locations whose time series are given in Fig. 2; see also Lambert and Boer, 2001, for discussion).
The comparisons above suggest that climate models may be informative about changes in climate, even while fail-ing to capture certain current characteristics.This is welldemonstrated for means (again, see Fig. 1), and the fact that AOGCMs capture trends in mean climate well suggests that their physics may be sufficiently realistic to provide a guide to trends in variability.We therefore seek a method of producing simulations of future climate that combines model output with data to incorporate both observational ground truth and model forecasts of trends.An appropriate method should simply reproduce current climate when models suggest no changes.When models do predict changes, the desired "data-driven simulation" should reproduce model changes in second-order moments (e.g., covariance) of climate but retain most non-Gaussian characteristics of data, rather than of model output, when changes in variability are relatively small.Our motivation in this work is to develop an empirically driven approach to simulating future climate that modifies existing observations in terms of means and secondorder moments (including covariances) based on changes in model simulations.
Many methods for combining observations with model output in climate projections have been developed for use in impacts studies, especially those involving hydrology and agriculture (see, e.g., Wood et al., 2004;Diaz-Nieto and Wilby, 2005;Eisner et al., 2012;Hawkins et al., 2013).In these cases, impacts models typically require inputs of temperature and precipitation at finer spatial resolution than is provided by AOGCMs, whose typical state-of-the-art resolution is on the order of 1 • (111 km or 69 miles).For this reason, approaches for simulating future climate by combining model output and data are often intertwined with methods for downscaling to higher spatial resolutions, and are described in the literature on statistical downscaling1 .We provide a  brief summary of existing approaches, along with what we consider to be the primary shortcomings of each approach.All approaches that combine observations and model output in simulating future climates correct in some way for model-observation discrepancies.One approach is a simple "bias correction" in which any offsets between current observed and modeled present-day climate are assumed to be systematic model errors.Model simulations of future climate are then "corrected" by adding the present-day bias (determined by comparing observations to a baseline run).Bias corrections can be made on annual mean temperatures or, more commonly, on monthly mean temperatures or annual harmonics, since models may not perfectly capture observed seasonal variation.One drawback of this approach is that all higher-order moments of the marginal and joint probability distributions (variability, skewness, stationarity, etc.) are provided by the future model output.As we have seen in Fig. 3, climate models may not adequately capture higherorder characteristics in the data.
A variant on this approach, typically termed "bias correction/spatial disaggregation" (BCSD), attempts to provide a better approximation of observed climate distributions by separately bias-correcting the different quantiles of model output (e.g., Wood et al., 2002Wood et al., , 2004)).This approach is also termed "quantile mapping" and involves computing a transfer function between model simulations of present-day climate and actual observations based on the ranked model output.The transfer function is then applied to AOGCM projections of future climate.This approach accommodates errors in higher-order moments of the model -in the most extreme case, the procedure results in a full transformation of the empirical cumulative distribution function (CDF) -but only corrects the marginal distributions of the model, and takes no account of differences in the covariance structure of the model output and the observed climate.Since human societies are sensitive to climate variation at different timescales (e.g., to changes in duration of droughts or rainfall that produces flooding), BCSD is not ideal for estimating the societal impacts of climate change.
While the previous two approaches are model-based, i.e., they quantify present-day model-observation discrepancy and apply it to future model output, the "change-factor" or "delta" method (see, e.g., Diaz-Nieto and Wilby, 2005;Eisner et al., 2012, and references therein) is observation-based: future climate projections are generated by modifying observations based on present-future differences in AOGCM simulations.Specifically, the delta method involves quantifying the difference (e.g., in annual or monthly means) between model output from a baseline run driven with present-day or preindustrial GHG concentrations and that from a "scewww.adv-stat-clim-meteorol-oceanogr.net/1/1/2015/Adv.Stat.Clim.Meteorol.Oceanogr., 1, 1-14, 2015 nario" run under future GHG concentrations 2 , then adding this difference onto some observation set.As a result, higherorder moments (in terms of the marginal and joint PDFs) will be derived from observations.Hawkins et al. (2013) showed that delta-method approaches may provide a better forecast of future climate than bias-correction approaches.Deltamethod approaches do not however generally involve representing changes in variability in future climate regimes.Recent advances have been developed to accommodate changing marginal variances (see, e.g., Ho et al., 2012;Hawkins et al., 2013); however, such approaches ignore joint dependence characteristics (e.g., covariances).
In this work, we adopt the observation-based approach of the delta method (modifying observations based on changes suggested by model output) but extend the method to account for possible changes in variability and temporal correlations.While recent work has extended the delta-method approach to accommodate some aspects of changing variability (Ho et al., 2012;Hawkins et al., 2013), these methods do not account for changes in third-order or higher moments of the marginal distribution, or in the covariance of the joint distribution.Changing covariance structures in particular are a critical component of simulating future climate for impacts assessment.
A delta-or change-factor approach that involves modifying covariance structures poses substantial challenges.The approach requires modifying a vector of random variables with a given joint dependence structure to produce a new vector of random variables with a different dependence structure.To achieve this goal, it helps to think about modifying quantities that are independent (or close to independent) under both present and future climates.In this regard, spectral-based approaches provide a natural framework.We propose an approach that modifies the discrete Fourier transform (DFT) of observations based on an estimated ratio of spectral densities of model output.Under a large class of stationary processes, the DFT is a transformation to approximate independence (Brillinger, 1981).This approach shares an important quality with the delta method that when the model suggests no changes (in either first-or second-order moment characteristics), the simulations equal the observations.
One caveat is that the procedure is designed to transform model simulations of an assumed equilibrium climate to another equilibrium climate while, during foreseeable human timescales, climate will continue to remain in a transient state.This approach does not directly address the important problem of simulating transient climate behavior in the covariance structure.However, it is likely that the method would remain an improvement over the delta method even in predicting future transient climate states, with certain exten-2 Because current climate is transient and changes as a result of increased GHG emissions are not fully realized yet, preindustrial GHG forcings may be a reasonable assumption in these problems.
sions related to nonstationary time series.We do not explore the issue in this paper, but point out a potential approach in Sect. 4.
In the remainder of the paper, Sect. 2 outlines the methodology, explaining how to estimate the ratio of spectral densities and use it to modify observations.We also explain an approach to account for a limited type of temporal nonstationarity in the data as brought about by differences in intraseasonal variability across seasons.Section 3 applies the method to generating simulations of daily mean temperature for a higher-CO 2 world, and Sect. 4 discusses results and future research needs.We provide supplemental materials that give further details, a numerical study, and information on how to access the code and data used to reproduce the analysis.

Methods
Our method produces data-driven simulations of future climate that combine observed climate with model predictions of changes to climate means, variability and temporal correlation.To do this we need to take account of changes in variability of model output over all temporal scales.
In the sections below, we first demonstrate the principle of our approach for an idealized situation: we assume an infinite length observational time series with known changes in the spectral process.We then develop the method for the more practical setting in which -the time series of both observations and model output are finite -we do not know the explicit form of the spectral process -we do not know the explicit form of changes to the spectral process -climate exhibits a strong seasonal cycle in both first and second-order moments.

Motivation
We demonstrate here that given an infinite length Gaussian time series representing present-day climate with a known spectral process and known future changes in the spectral process, we can modify the continuous Fourier transform separately at each frequency to produce output that has the correct joint distribution for the future process.Let {Z 0,t ; t = 0, ±1, ±2, . ..} represent a time series of an observable process of interest.Furthermore, suppose {Z 0,t } is a stationary Gaussian process with E(Z 0,t ) = 0 and covariance function Let {Z 1,t ; t = 0, ±1, ±2, . ..} represent the future process that we wish to simulate.Suppose that {Z 1,t } is also a stationary Gaussian process with We are interested in modifying {Z 0,t } in order to generate a random process that is equal in (joint) distribution to {Z 1,t }.The temporal correlations in {Z 0,t } makes this nontrivial.However, the orthogonal nature of the spectral representation makes it the natural domain in which to modify random quantities.For example, writing ı for √ −1, Z 0,t has the representation Z 0,t = 0.5 −0.5 exp(2π ıωt)d Ẑ0 (ω) where Ẑ0 (ω) is a complex-valued Gaussian random measure with mean of 0 and for disjoint sets ([c, d])) = 0, and x represents the complex conjugate of x.That is, Ẑ0 (ω) (also referred to as the spectral process associated with {Z 0,t }) is an orthogonal Gaussian measure that can be modified separately at each frequency to produce a process with a different covariance structure.The spectral distribution associated with {Z 0,t }, G 0 (ω), is the positive finite measure given by E( d Ẑ0 (ω) 2 ) = dG 0 (ω).Assuming absolute summability of the covariance, i.e., spectral distribution is absolutely continuous: dG 0 (ω) = g 0 (ω)dω and g 0 is called the spectral density for {Z 0,t }.
The spectral density can be obtained from the covariance summable, we can similarly consider the second-order characteristics of {Z 1,t } based on its spectral density The spectral densities g 0 (ω) and g 1 (ω) provide information regarding the covariance structure of {Z 0,t } and {Z 1,t }, respectively, for frequencies, ω ∈ (−0.5, 0.5].Then, given a spectral density ratio ρ g (ω) = g 1 (ω)/g 0 (ω) we can modify the spectral process associated with Z 0,t to generate which is a stationary, Gaussian process with E(Z 1,t ) = 0 and covariance γ 1 (h).In this way, we derive the future, unobservable process in terms of the present process, modified by the ratio of their spectral densities.If g 1 (ω) = g 0 (ω), for all ω ∈ (−0.5, 0.5], then Z 1,t = Z 0,t , for all t (because in this special case the procedure reduces to taking the DFT and then the inverse DFT of the observations).In particular, the temporal covariance structure of the simulations equals that of the observations.

Outline of approach
When working with real time series of climate observations and model output, the spectral densities in the past and future, g 0 (ω) and g 1 (ω), are not known.While g 0 (ω) can be estimated from data, clearly we cannot provide an observationbased estimate of g 1 (ω).A central question then becomes how to best represent the spectral ratio ρ g (ω).Let f 1 (ω) represent the future spectral density associated with the computer model output.For AOGCMs, f 1 (ω) may differ substantially from g 1 (ω).However, given the model's suggested covariance structure under a baseline period, represented by f 0 (ω), the estimated change in covariance structure may be a reasonable approximation for the real changes in the covariance structure, especially if those changes are relatively small.We therefore do not assume that model output has the correct covariance structure for a given GHG scenario, but assume that the computer model provides a reasonable approximation to the changes in the spectral density across all frequencies (i.e., ρ g (ω) = ρ f (ω) for all ω ∈ (−0.5, 0.5], where Carrying out the simulation on real data then requires the following steps, starting with {Z 0,t } (observations), {Y

Reverse preprocessing to produce simulations Z
In the following subsections we describe in detail these primary steps: estimating the spectral ratio and modifying the discrete Fourier transform; removing the seasonal cycle; and modulating the deseasonalized time series.

Spectral-based conditional simulation
Let {Z 0,t ; t = 0, . .., T − 1} represent the observations of the process of interest, observed at regular time points.For now, assume that the process is stationary with E(Z 0,t ) = 0. We discuss how we account for any trend in Sects.2.5 and 2.6.
In the previous section T = ∞ whereas, in practice, our observations are observed discretely over a finite period and the spectral process associated with Z 0,t is unknown.First we approximate the true spectral process by using the discrete Fourier transform of the obser- and also Cov( Ẑc 0,j , Ẑs s,k ) = T , for all j, k. (See Shumway and Stoffer, 2011, for details.) Here, T is a generic remainder that varies with j and k and can be shown to obey a bound | T | ≤ β/T .As a result, our methodology is modifying nearly independent quantities in order to produce simulations with a different covariance structure than the observations.Let ρf (ω k ) be an estimate of the ratio of the spectral densities at ω k .Then, the simulations (not accounting for changes in mean) under a given scenario can be represented as (1) so, when ρf This suggests the following covariance structure for {Z 1,t } for a given estimate ρf (ω k ): A brief numerical study in Sect.S3 illustrates the efficacy of this approach even for fairly small T when ρ g = ρ f is known.In the following section, we provide the details of a penalized likelihood approach to estimate ρ f (ω k ).Finally, although we have motivated this methodology in terms of Gaussian processes, the resulting simulation of Z 1,t in Eq. (1) will tend to retain any non-Gaussian characteristics of Z 0,t , at least if ρf (ω k ) is nearly constant.

Estimation of the ratio of spectral densities (ρ f ω k )
We propose a penalized likelihood approach for estimation of ρ f (ω k ), similar to the approach given in Pawitan and O'Sullivan (1994) for the estimation of one spectral density.Let f j,k = f j (ω k ), θ j,k = log(f j,k ), j = 0, 1, k = 1, . .., K, and θ j = (θ j 1 , . .., θ j K ) .Then, a penalized likelihood can be generally written as where L 0 (θ 0 ) and L 1 (θ 1 ) represent the Whittle likelihood for j = 0, 1, respectively.So, L 0 (θ 0 ) and L 1 (θ 1 ) provide an objective function that determines the fit to the data, J (θ 0 , θ 1 ) is a function that penalizes lack of smoothness, and δ is a smoothness parameter.
Let Ŷi,j,k = T −1/2 T −1 t=0 Y i,j,t exp(−2π ıω k t) represent the DFT of the ith realization of the model output at frequency ω k (for either the model baseline or scenario run).Note that when Y i,0,t and Y i,1,t follow stationary, Gaussian distributions, the periodograms I i,j,k = Ŷi,j,k 2 , for j = 0, 1, follow (asymptotically) independent exponential distributions such that E I i,j,k /f j,k → 1 as T → ∞.As a result, the Whittle negative-log-likelihood approximation L(θ j ; I i,j ) = T /2 k=−T /2+1 θ j,k + I i,j,k exp(−θ j,k ) is a reasonable approximation for the likelihood in the objective function (Whittle, 1954).In the event that we have multiple realizations, we can take the average periodogram I j,k = n j i=1 I i,j,k which follows asymptotically a gamma distribution with E(I j,k ) = f j,k as before but with Var(I j,k ) = f j,k /n j (as opposed to Var(I j,k ) = f j,k ).
We further linearize the log likelihood and carry out estimation using an iterative, weighted least squares approach (McCullagh and Nelder, 1989) Thus, this framework can handle the situation in which there is a different number of independent realizations for the baseline and scenario runs.For given initial conditions, these computations are iterated until convergence.

Penalty function
Although penalties could be placed on the individual spectral densities themselves, for our analysis we only need an estimate of the ratio; therefore, we place the penalty on the log ratio of the spectral densities θ 1 − θ 0 so that J (θ 0 , θ 1 ) can be written as J (θ 1 − θ 0 ).Because we expect the ratio of spectral densities to be smoother than the individual spectral densities themselves, it makes sense to place the penalty on this ratio, enabling us to obtain a low-variance estimate of the ratio while increasing the bias less than we would by smoothing each spectral density individually.Our penalty function is then placed on the th derivative of λ(ω) = θ 1 (ω) − θ 0 (ω): J (λ) = (2π ) −2 0.5 −0.5 λ ( ) (ω) 2 dω.
Using Parseval's identity, this can be written as , where * k is the kth Fourier coefficient of λ(ω), * k = 0.5 −0.5 λ(ω) exp(−2π ıkω)dω.We then approximate the penalty function by where k is the discrete Fourier coefficient of λ, k = The objective function that we minimize can then be written as where, for a given smoothing parameter δ, we can iterate back and forth between estimates of θ 0,k and λ k until convergence.The ratio of spectral densities can be estimated using the algorithm provided in Sect.S1.We do not develop an automated method for choosing the smoothing parameter δ in this paper.In a situation in which multiple realizations of a climate scenario exist, it may be desirable to choose δ based on a cross-validation study.Here, we chose δ = e −7 ≈ 9.12×10 −4 , which appears to give good visual results.Using the formula for effective degrees of freedom given by Pawitan and O'Sullivan (1994) yields an approximate bandwidth for this smoother of 0.078 day −1 , which is quite broad considering that we are defining the spectral density on (−0.5, 0.5] day −1 .We believe this degree of smoothing is acceptable given that the estimated log-spectral ratios are quite flat.(As mentioned previously, one advantage of smoothing on the ratio of spectral densities is that the ratio is flatter than are the individual spectral densities.)However, we do see some evidence that the ratio is less flat at the lowest (belowannual) frequencies.For studies of interannual variability, there could be some advantage in using a penalty function that allows for more flexibility in λ near the origin by defining J (λ) = (2π ) −2 0.5 −0.5 η(ω){λ ( ) (ω)} 2 dω for some positive, even function η that takes on smaller values near 0. Such a technique could resolve different changes at different interannual frequencies.

Seasonal cycle and long-term trend
The previous section assumed that the process of interest was a stationary process with constant mean.Daily mean temperature however involves a strong seasonal component.So, before estimating the spectral ratio and modifying the DFT of the observations, we remove the seasonal cycle in the observations and AOGCM output.The empirical mean of the observations and present-future difference in the AOGCM output are then added back on at the end of the algorithm.This part of our approach is analogous to the delta method and in fact reduces to the delta method when the present and future spectral densities are equal.
As mentioned previously, the delta method uses model output for changes in first-order characteristics (e.g., overall mean and seasonal cycle) estimated from the model output.This method typically involves adding the difference in the overall mean (usually including the seasonal cycle) of the base and scenario time slices for the AOGCM to the observations.Let μ0,t and μ1,t represent monthly means or annual harmonics, i.e., μj,t = μj + K k=1 Rj,k cos(2π ω k t + φj,k ) for j = z, 0, 1, and ω k = k/365.25.The parameters μz , μ0 , and μ1 are the estimated long-term averages for the observations, base, and scenario periods, respectively, and Rz,k , R0,k , and R1,k are the estimated amplitudes at ω k for the observations, base, and scenario periods, respectively.Lastly, φz,k , φ0,k , and φ1,k are the estimated phase shifts for ω k .All parameters are estimated using least squares.Seasonal demodulation (Sect.2.6) is performed on Z 0,t = Z 0,t − μz,t , Y 0,t = Y 0,t − μ0,t , and Y 1,t = Y 1,t − μ1,t , in order to account for seasonal difference in second-order moments.In our example, our AOGCMs have been run far past the point of CO 2 stabilization and, therefore, can be considered to be nearly in an equilibrated state.However, there is evidence of a long-term trend in temperature in the observations {Z 0,t }.We remove this long-term trend from the observations using a simple linear regression of the observations against the logarithm of CO 2 3 .

Accounting for seasonal nonstationarity
Thus far we have assumed that the deseasonalized observations and model output, Z 0,t = Z 0,t − μz,t , Y 0,t = Y 0,t − μ0,t , and Y 1,t = Y 1,t − μ1,t are (temporally) stationary.However, this need not be the case and in applications involving daily mean temperature it likely is not the case.Figure 4 shows the log-averaged periodograms by season for the Illinois pixel for the base and scenario period, as well as for the observations (similar plots are provided for the Southern Ocean 3 The trend in the observations may be affected by volcanoes (e.g., Pinatubo), which produce a temporary reduction in GMT.The fact that these trends are not removed implicitly assumes that intermittent volcanic eruptions would continue in the future.Another potential concern is that the aerosol forcings that affect observed climate will not continue to evolve indefinitely as they have in the past.
and Gulf of Guinea in Figs.S1 and S2 in the Supplement).Clearly, the seasonal spectral density functions for the base period are different for the different seasons, with the winter months showing the greatest variability and the summer months the lowest variability across all frequencies.Note that in the case of the scenario period, variability across frequencies in the winter has decreased, and is now roughly the same as in spring and fall, but the summer variability is still roughly the same, and is lower across most frequencies.Thus, the assumption of temporal stationarity is not reasonable and, furthermore, the form of the nonstationarity is somewhat different for the base and scenario periods.However, the log periodograms for the different seasons are nearly parallel for both periods, suggesting that it may be reasonable to treat the processes as uniformly modulated (Priestley, 1988).
Following Priestley, we consider , and Y 0,t = D 0,t Y * 0,t , after deseasonalizing, where Z * t , Y * 1,t , and Y * 0,t are stationary processes (corresponding to the observations, model output under scenario period, and model output under base period, respectively).Then, D z,t , D 1,t , and D 0,t are modulation constants to be estimated.Thus, if we can find suitable values for the modulation constants, then we can perform the spectral den-sity estimation on Y * 1,t and Y * 0,t , in order to modify the DFT of Z * 0,t , and then multiply by the constants D 1,t as a last step to account for the nonstationarity across seasons.Our approach to estimating the modulation constants is provided in Sect.S2. Figure 4 and Figs.S1 and S2 show the logaveraged periodograms for the modified process.The logperiodograms are much closer together than they were originally, suggesting that this approach accounts for most of the seasonal nonstationarity.

Application
In this section, we continue to illustrate our methodology using NCEP Climate Forecast System Reanalysis observations (Saha et al., 2010) and the CCSM3 output described in Appendix A. Because the observations and model output used are of different lengths (32 years and 100 years, respectively) the Fourier frequencies will be different.As a result, after estimating the ratio of spectral densities for the model output, we do a simple linear interpolation on the log-spectral ratio of the model output to the Fourier frequencies of the observations.
Although we are only modifying the temporal covariance structure, we can produce maps that show how variability is changing at different locations and different frequencies (e.g., see Fig. 5).In general, at most locations and at most frequencies, variability is decreasing in the CCSM3 output for this particular scenario.Variability increases occur only in a few regions.Increases in lower-frequency (periods close to 3.2 years) variability appear primarily on land at low latitudes.Increases in higher-frequency (periods of roughly 2 days) variability occur primarily at low latitudes over water near coastlines.
Variability clearly changes differently at different locations (Figs. 5, 6) and, furthermore, variability changes at a given location can differ with frequency (Fig. 6, top panel).In Illinois and the Gulf of Guinea, there is a modest decrease in low-frequency variability.At high frequencies, there is a slightly greater suppression of variation in Illinois, whereas in the Gulf of Guinea high-frequency variation is actually larger for the future scenario than the present.The decreasing variability at high frequencies in Illinois may be consistent with suggested changes in the polar jet stream that impacts weather at the middle latitudes.For the Gulf of Guinea, the slight suppression of low-frequency variation and the amplification of high-frequency variation may suggest fundamentally changing weather patterns at this location.These results show that the manner in which variability changes is nontrivial and is dependent on the temporal scale.As a result, an approach that considers changes across a variety of timescales is necessary (as opposed to a simple rescaling of the observations based on changes in model output).
In contrast to those locations, however, are pixels such as the Southern Ocean, where the change in variability re- The low-frequency results are the estimated log ratios at 1168 days and the high-frequency results at 2 days; however, due to the large degree of smoothing, it is best to think of them as lowand high-frequency behavior.Both long-(left column) and shortterm (right column) variability decreases in nearly all locations.Remainder of rows: estimated log-spectral densities at these frequencies for reanalysis (second row), model baseline period (third row) and model scenario period (bottom row), using the demodulated and deseasonalized time series.The pattern of enhanced variability over land vs. ocean and high vs. low latitudes is as expected.
mains relatively constant across all frequencies (with approximately a 60 % decrease in overall variability).For locations that exhibit this type of change in the spectral ratio, a simple scaling of the observations may be acceptable.However, Fig. 6 indicates that in all three locations, the acrossfrequency variation of the spectral density is greater than the across-frequency variation of the spectral ratio, supporting our claim that the spectral ratio is smoother than the spectral densities themselves.
All three locations used as examples show evidence of a mean shift (see Fig. ) of the estimated spectral density ratios in the Southern Ocean (blue), Illinois (green), and Gulf of Guinea (red).(Bottom) Logarithm of the estimated spectral densities of the reanalysis data (solid line), base period (dashed line), and scenario period (dashed and dotted line).The spectral density estimation was performed on the deseasonalized and demodulated time series.
is small relative to temperature variability.In the Southern Ocean location, the mean shift in local Winter (JJA, June-August) is nearly 10 • , which is likely due to the loss of sea ice in the future scenario.(Ice cover allows for lower temperatures than are possible over open ocean).All locations show physically reasonable characteristics in variability and in changes in variability: variability is stronger in the the nonequatorial locations (Southern Ocean and Illinois) than near the Equator and stronger in winter than in summer, and variability reductions are also greater in winter.
An important aspect of our approach is that it does not significantly alter the non-Gaussian aspects (e.g., tail behavior) of observed climate.In fact, in our method, when the model does not show changes in mean or covariance, the simulations are simply the observations and, thus, non-Gaussian features of the data are retained exactly.When the observations are not significantly changed, the non-Gaussian features of the data are largely retained.For instance, JJA in the Gulf of Guinea shows a marginal distribution that is positively skewed.In this case, the simulation shows a slight decrease in marginal variability, as well as an increase in mean temperatures, but we maintain the positive skewness of the observations (see Fig. 8).We consider this to be a strength of our approach: in the event that there are non-Gaussian features of the data, the simulations will retain these features, at least when the change in the dependence structure is relatively small.
Preserving the shapes of distribution of the observations (e.g., skewness, kurtosis) would be a problem if the actual shapes of distributions changed from present to future.For locations in or near bodies of water, changes in temperature means can alter climate variability distributions because those distributions are sensitive to the freezing point of wa- ter.As long as water is liquid, temperature variability is constrained because air temperatures cannot drop significantly below freezing.This property is evident in the marginal distributions for both observations and model output in Fig. 3, which show upward bumps in the marginal densities at nearfreezing temperature, i.e., a clustering of temperatures near freezing.When open water freezes, however, air temperatures can become very cold.The left-skewed tail in the Southern Ocean local winter in reanalysis data (Fig. 3, bottom row JJA) indicates that the location is sea-ice covered for part of the season.This freezing point effect introduces further problems when model output is biased relative to observations.The example of the Southern Ocean location in winter shows this clearly.The wintertime base period CCSM3 temperatures show a distribution characteristic of complete sea ice cover rather than the partial cover implied by observations, with a strong cold bias (a mode over 10 • below the freezing point) and a wide spread.Because sea ice is lost in the warmer future scenario, the temperature change between model base and future periods is very large.The modeled changes in both mean temperature and temperature variability are therefore physically inappropriate when applied to observations and produce unphysically high wintertime temperatures and low wintertime variability in resulting data-driven simulations (Fig. 7).We note that this limitation applies not just to our approach but to statistical downscaling methods in general.Systematic offsets in the mean between climate model and data make it difficult to adequately simulate future climate, so that simulations are inherently limited by the skill of the AOGCM being used.for the pixels at Illinois (top row), Gulf of Guinea (middle row), Southern Ocean (bottom row) for reanalysis (blue line) and the simulations (red line).The simulations display marginal tail behavior similar to the reanalysis observations.

Discussion
Detailed characterization of the nature in which climate is changing (mean shifting, tail behavior, spatial and temporal covariance structures) is still a relatively open area of inquiry.One of the best ways of studying how future climate might change is by first investigating the nature in which the statistical properties of the output of AOGCMs change from the present to possible future scenarios.We have provided a method of quantifying how temporal covariance is changing in these AOGCMs at different temporal scales.Our results show that variability is changing differently at different locations.At a given location, the changes in variability may be different (in both magnitude, direction of the change, or both) across different frequencies.We used this estimate of changing variability to produce simulations that modify the temporal covariance structure of the observations.In this way, we extend the delta method to be able to account for changes in the mean and covariance structure.
Our for producing simulations relies on modifying the discrete Fourier transform of the observations and, as such, the length of the simulations in this manner is currently limited by the available data.However, it is possible that by recycling old observations, one could generate simulations of longer length.Another possibility is to modify the observations by phase scrambling (Davison and Hinkley, 1997) and then appending these newer pseudo-observations to the true observations.
We point out that we have not accounted for any changes in spatial and spatiotemporal covariance structures.Accounting for changes in spatial covariance is complicated by the nonstationarity present in the observations (due to geogra-phy, land-ocean contrasts, etc.) and remains the subject of future research.However, we do note that, due to the use of the observations, we are mimicking any spatial structure in the present climate regime.
Next, while we have provided a method for producing simulations of daily mean temperature, most impacts assessments also rely on simulations of daily precipitation.The methodology presented here is not fit to handle the highly non-Gaussian, nonlinear nature of daily precipitation directly; however, a popular approach in the statistics literature is to model precipitation using a latent Gaussian process (see Allcroft and Glasbey, 2003, for an example).The approach presented could be applied to such a latent process.Latent processes might be further extended to consider joint bias correction and downscaling of temperature and precipitation.
Perhaps most importantly, the methodology presented here is based on the assumption of stationarity in the model output and the data.While we did incorporate the concept of a uniformly modulated process to deal with seasonal nonstationarity, this methodology is still limited to simulating equilibrium climate.Because for the foreseeable future our climate will be in a transient state, we must consider ways of extending this methodology to be able to simulate transient climate.We point out that there is the potential for this methodology to be extended by considering an evolutionary spectral approach (Priestley, 1988).
Lastly, our methodology is limited to generating simulations for those GHG scenarios for which an AOGCM has been run.We cannot produce simulations for arbitrary GHG emissions scenarios without first running the AOGCM to obtain the necessary output.However, we note the potential to consider "emulating" higher-order characteristics in the general circulation models in order to generate simulations for arbitrary emissions scenarios.For transient climates, it may be possible to relate changes in the covariance structure to the past trajectory of CO 2 .
We believe that our approach provides a general framework for high-resolution future climate simulation.Two critical features of our approach are that (1) it is observationdriven, using the model output to suggest how to modify the existing observations, and (2) it considers changes in both mean and covariance structures; and this modification of covariance structure, by taking place in the frequency domain, involves modifying quantities that are at least approximately independent.We anticipate many opportunities to extend our framework to generate more realistic simulations for use in impacts assessment, and suggest that any extensions should seek to maintain these features when feasible.Society's obvious need for better impacts assessment begins with a better understanding of how climate will change in the future.Our conceptual approach provides a valuable framework for quantifying climate change and simulating future climate in order to meet that need.

Figure 1 .
Figure 1.Comparison of modeled and observed global mean temperatures.(a) CO 2 concentrations used in "baseline" and "scenario" runs with the CCSM3 model.(See Appendix A for description of experiments and observations.)Figures here truncate output after less than 600 years but the scenario run extends for 6000 years.(b) Corresponding annual model GMT ( • C) for the two runs and observed GMT from the Global Historical Climatological Network.Model output reproduces trends in global temperature well but with a systematic offset from observations.To better show the similarity in trend we also plot the observational record minus a 2 • C offset.

Figure 2 .
Figure 2. (Left) Three locations (individual model pixels) used as examples throughout the manuscript, chosen to represent different combinations of seasonality, variability, and expected future changes: Illinois, mid-continental with a strong seasonal component (green, 38.97 • N, 90 • W); Gulf of Guinea, near-equatorial with little seasonal cycle (red, 1.86 • S, 0 • E); and Southern Ocean, which has strong projected changes in both mean temperature and in variability (blue, 61.2 • S, 33.8 • E).Annual standard deviation of daily temperatures σ and projected temperature change (scenario-baseline) are Illinois: σ = 10.81,= 3.87; Gulf of Guinea: σ = 1.97, = 2.43; Southern Ocean: σ = 4.67, = 8.10.(Right) Time series of the 3 years of daily temperature ( • C) from the NCEP-DOE (National Centers for Environmental Predictions -Department of Energy) Climate Forecast System Reanalysis at those locations.(See Appendix A for a description of the observational data set.)

Figure 3 .
Figure3.Marginal densities (by season) of daily mean temperature ( • C) for the pixels in Illinois (top row), the Gulf of Guinea (middle row), and the Southern Ocean (bottom row) for reanalysis data (solid blue line), baseline model output (dashed blue line), and scenario model output (dashed red line).The model output does not replicate the marginal distributions of the reanalysis observations.Furthermore, the marginal distributions in the model output change from the baseline to scenario periods.

Figure 4 .
Figure 4. (Top) Log (base 10) of averaged periodograms for the Illinois location, by season, for the reanalysis (left), model baseline period (middle), and model scenario period (right).Note strongest variability in winter, weakest in summer.(Middle) Identical to top but now for the demodulated time series.Seasonal differences in variability are effectively removed, suggesting we can treat these time series as stationary in time.(Bottom) Modulation constants used for the reanalysis (left), model baseline period (middle) and model scenario period (right), showing smallest values in summer, as expected.See Figs.S1 and S2 for similar plots for other locations used as examples; results are similar.

Figure 5 .
Figure 5. (Top) Estimated log (base 10) ratio of spectral densities for model scenario vs. baseline at low and high frequencies.The low-frequency results are the estimated log ratios at 1168 days and the high-frequency results at 2 days; however, due to the large degree of smoothing, it is best to think of them as lowand high-frequency behavior.Both long-(left column) and shortterm (right column) variability decreases in nearly all locations.Remainder of rows: estimated log-spectral densities at these frequencies for reanalysis (second row), model baseline period (third row) and model scenario period (bottom row), using the demodulated and deseasonalized time series.The pattern of enhanced variability over land vs. ocean and high vs. low latitudes is as expected.

Figure 6 .
Figure6.(Top) Logarithm (base 10) of the estimated spectral density ratios in the Southern Ocean (blue), Illinois (green), and Gulf of Guinea (red).(Bottom) Logarithm of the estimated spectral densities of the reanalysis data (solid line), base period (dashed line), and scenario period (dashed and dotted line).The spectral density estimation was performed on the deseasonalized and demodulated time series.

Figure 7 .
Figure 7. Time series of daily mean temperature ( • C) for 3 years of reanalysis (blue) and simulations (red) at Illinois (top row), Gulf of Guinea (middle row), and Southern Ocean (bottom row).

Figure 8 .
Figure 8. Marginal densities (by season) of daily temperature ( • C) for the pixels at Illinois (top row), Gulf of Guinea (middle row),Southern Ocean (bottom row) for reanalysis (blue line) and the simulations (red line).The simulations display marginal tail behavior similar to the reanalysis observations.
0,t } (model base period time series), and {Y 1,t } (model scenario period time series): )/D 1,t , which have mean of 0 and are stationary.See Sect.2.5 for details on the estimation of the seasonal cycle (i.e., μz,t , μ0,t , and μ1,t ) and see Sect.2.6 and Sect.S2 in the Supplement for details on the estimation of seasonal variation (i.e., D z,t , D 0,t , and D 1,t ).