A method to encapsulate model structural uncertainty in ensemble projections of future climate : EPIC v 1 . 0

A method, based on climate pattern scaling, has been developed to expand a small number of projections of fields of a selected climate variable (X) into an ensemble that encapsulates a wide range of indicative model structural uncertainties. The method described in this paper is referred to as the Ensemble Projections Incorporating Climate model uncertainty (EPIC) method. Each ensemble member is constructed by adding contributions from (1) a climatology derived from observations that represents the time-invariant part of the signal; (2) a contribution from forced changes in X, where those changes can be statistically related to changes in global mean surface temperature (Tglobal); and (3) a contribution from unforced variability that is generated by a stochastic weather generator. The patterns of unforced variability are also allowed to respond to changes in Tglobal. The statistical relationships between changes in X (and its patterns of variability) and Tglobal are obtained in a “training” phase. Then, in an “implementation” phase, 190 simulations of Tglobal are generated using a simple climate model tuned to emulate 19 different global climate models (GCMs) and 10 different carbon cycle models. Using the generated Tglobal time series and the correlation between the forced changes in X and Tglobal, obtained in the “training” phase, the forced change in theX field can be generated many times using Monte Carlo analysis. A stochastic weather generator is used to generate realistic representations of weather which include spatial coherence. Because GCMs and regional climate models (RCMs) are less likely to correctly represent unforced variability compared to observations, the stochastic weather generator takes as input measures of variability derived from observations, but also responds to forced changes in climate in a way that is consistent with the RCM projections. This approach to generating a large ensemble of projections is many orders of magnitude more computationally efficient than running multiple GCM or RCM simulations. Such a large ensemble of projections permits a description of a probability density function (PDF) of future climate states rather than a small number of individual story lines within that PDF, which may not be representative of the PDF as a whole; the EPIC method largely corrects for such potential sampling biases. The method is useful for providing projections of changes in climate to users wishing to investigate the impacts and implications of climate change in a probabilistic way. A web-based tool, using the EPIC method to provide probabilistic projections of changes in daily maximum and minimum temperatures for New Zealand, has been developed and is described in this paper.

T global time series and the correlation between the forced changes in X and T global , obtained in the 'training' phase, the forced change in the X field can be generated many times using Monte Carlo analysis. A stochastic weather generator model is used to generate realistic representations of weather which include spatial coherence. Because GCMs and Regional Climate Models (RCMs) are less likely to correctly represent unforced variability compared to observations, the stochastic weather generator takes as input measures of variability derived from observations, but also responds to forced changes in climate in a way that is 15 consistent with the RCM projections. This approach to generating a large ensemble of projections is many orders of magnitude more computationally efficient than running multiple GCM or RCM simulations. Such a large ensemble of projections permits a description of a Probability Density Function (PDF) of future climate states rather than a small number of individual story lines within that PDF which may not be representative of the PDF as a whole; the EPIC method corrects for such potential sampling biases. The method is useful for providing projections of changes in climate to users wishing to investigate the impacts 20 and implications of climate change in a probabilistic way. A web-based tool, using the EPIC method to provide probabilistic projections of changes in daily maximum and minimum temperatures for New Zealand, has been developed and is described in this paper.

Introduction
While future changes in climate will follow a single trajectory, it is highly unlikely that any single climate model projection will simulate that trajectory. The use of a single model projection is therefore insufficient for assessing the potential future state of the climate. Rather, what is required is a large (e.g. 10,000 member) ensemble of projections that provides a probabilistic portrayal of how the climate is expected to evolve. Clustering of trajectories within that probabilistic envelope then shows where 5 any single trajectory has a higher likelihood of occurring. Probabilistic simulations of future climate, presented as Probability Density Functions (PDFs), give decision-makers a much clearer picture of likelihoods of certain future climate states compared to a single projection, or a small set of projections (Watterson et al., 2008). That said, if decision-makers are presented with PDFs obtained from the same family of models, these may be biased by the assumptions and limitations inherent in a single family of models that do not explore the possible trajectories seen in other model families. PDFs of future climate that consider 10 all sources of uncertainty, including uncertainty resulting from structural differences in the underlying models, provide the information needed for quantitative risk assessments, since the likelihood of any particular trajectory can be estimated.
Exploring expected changes in extreme weather events also requires probabilistic simulations of future climate. While climate change may result in a small shift in the mean and/or standard deviation of a PDF of a selected climate variable, the tails of the distribution, which represent extreme weather events, can exhibit fractionally much larger changes (see Figure 1.8 in 15 Solomon et al. (2007)). It is especially important that extreme events, which by their nature are unusual, are captured in an ensemble of projections.
Resolving changes in the frequency of regional-scale extreme weather events requires large ensembles of projections of high spatial and temporal resolution. Generating such ensembles using models which simulate all important physical processes, such as Global Climate Models (GCMs) or Regional Climate Models (RCMs), is currently computationally prohibitive. The 20 ideas underlying climate pattern-scaling suggest a means of overcoming this hurdle and form the basis for the newly developed Ensemble Projections Incorporating Climate model uncertainty (EPIC) method described here. First, a robust statistical relationship is derived between the climate variable of interest (X) and some associated readily generated predictor. In climate pattern-scaling, this predictor is typically the global mean surface temperature (T global ). If observations are being used to establish this relationship, then observed values of X and T global would be used. If GCM or RCM output is used to establish the 25 relationship, then X and T global should come from the same model simulation.
Once the relationship between X and T global has been determined, then, given multiple versions of T global , multiple time series of X can be generated based on that relationship. This methodology assumes that many versions of T global can be simulated in such a way as to capture the inherent variability resulting from structural uncertainties in GCMs and carbon cycle models in a computationally efficient way, e.g., through the use of a simple climate model (SCM). If the large ensemble of 30 T global time series spans the range of model structural uncertainties, then the resultant ensemble of generated X time series will reflect that spread in uncertainties, e.g., as done in Reisinger et al. (2010).
2 Models and Data Sources 2.1 Regional Climate Model An RCM simulation, or a number of RCM simulations, are used to provide the time series used to train EPIC i.e. to quantitatively establish the relationship between the change in annual mean global mean surface temperature and the change in the climate variable of interest and its variability. The RCM simulations used in this study were performed using the Hadley Centre 5 regional climate model HadRM3-PRECIS (Jones et al., 2004) that has been modified to be used for New Zealand (Bhaskaran et al., 1999(Bhaskaran et al., , 2002Drost et al., 2007) Drost et al., 2007). The spatial resolution necessitates a computational time step of 3 minutes. The model orography and vegetation data sets were updated from those used by Drost et al. (2007) to the high resolution surface orography data set used in NIWA's operational forecast model (Ackerley et al., 2012); differences in the vegetation fields are small. The first year of model simulation (the spin-up) is excluded from the analysis as this is used to achieve quasi-equilibrium conditions of the land surface and the overlying atmosphere.  (Pope et al., 2000). Processes related to clouds, radiation, the boundary layer, diffusion, gravity wave drag, advection, precipitation and the sulfur cycle are all parameterized in HadAM3p. Additional details regarding HadAM3p are available in Gordon et al. (2000), Pope and Stratton (2002), Pope et al. (2000), and Gregory et al. (1994).
Most GCM and RCM simulations display biases when compared to observations. The RCM simulations used in this study 25 were partially bias-corrected by bias-correcting the sea surface temperatures that are used as lower boundary conditions for the HadAM3p simulations, which then provided the lateral boundary conditions for the RCM simulations (personal communication A. Sood, 2016).

Simple Climate Model
In this study, MAGICC (Model for Assessment of Greenhouse-gas Induced Climate Change; ; Mein-30 shausen et al. (2011)) is the simple climate model (SCM) used to generate an ensemble of T global time series. MAGICC is a reduced complexity climate model with an upwelling diffusive ocean and is coupled to a simple carbon cycle model that includes carbon dioxide (CO 2 ) fertilization and temperature feedback parameterisations of the terrestrial biosphere and oceanic uptake. MAGICC can be tuned to emulate the behaviour of 19 different AOGCMs (Meehl et al., 2007) and 10 carbon cycle models (Friedlingstein et al., 2006). The resultant 190 different 'tunings' for MAGICC can be used to generate 190 equally probable T global time series that provide some representation of the spread in T global resulting from structural uncertainties in AOGCMs and the carbon cycle models used in C4MIP (Coupled Carbon Cycle Climate Model Intercomparison Project).
When used as predictors for changes in local climate variables, and using the prior established quantitative relationship between 5 T global and the local climate variable, these 190 T global time series can be used to generate 190 equally probable time series of the local climate variable.

Virtual Climate Station Network
While the RCM simulations have been partially bias corrected, we recognise that some biases may remain. Therefore, we build our projections off an observational data set, so that, in the absence of any forced changes in climate, the projections default 10 to observations (this is described in greater detail below). Observationally-based time series are obtained from the so-called Virtual Climate Station Network (VCSN). The VCSN data set for the New Zealand land surface is constructed on a regular 0.05 • ×0.05 • grid from spatially inhomogeneous and temporally discontinuous quality controlled weather station data (Tait et al., 2005). The values estimated on the 0.05 • ×0.05 • grid are based on thin plate smoothing spline interpolation using a spatial interpolation model as described in Tait, (2008).

Methodology
For a given geographic location, each ensemble member, covering the period 1960 to 2100, is constructed from contributions including: 1. a climatology derived from observations that represents the time invariant part of the signal, 2. a contribution from long-term forced changes in the magnitude of the variable of interest where those changes scale with 20 changes in anomalies in global mean surface temperature (T global ), and 3. a contribution from weather generated by a stochastic weather generator that incorporates both forced and unforced variability.
The construction of each of these signals is described in greater detail below with a high level overview of how these contributions are related shown in Figure 1. The methodology described below pertains to a selected single greenhouse gas (GHG) 25 emissions scenario and the daily maps of the climate variable of interest (X; here daily maximum (T max ) and daily minimum (T min ) surface temperatures) are obtained from one or more RCM simulations. The results presented in this study were obtained from a 1900 member ensemble (10 Monte Carlo resampled alpha values were choosen for each of the 190 T g lobal time series from MAGICC), that was generated over the period 1960 to 2100. The T max and T min fields were obtained from five RCM simulations driven by the Representative Concentration Pathway (RCP) 8.5 GHG emissions scenario. RCP 8.5 was 30 choosen to be presented due to the high climate signal to noise ratio which results in more robust regression results, but the methodology is valid any choosen GHG emission scenarios. All anomalies were calculated with respect to the period 2000 to 2010.

The climatology
At each 0.05 • by 0.05 • (approximately 5km) grid point, a mean annual cycle is calculated from daily observational data from 5 2001 to 2010. For this study, these observational data were obtained from VCSN (Section 2.3). Since the 10 year baseline period is rather short, a climatology derived by calculating calendar day means would still contain some weather-induced noise. Therefore, a regression model which includes an offset basis function, expanded using two Fourier series expansions (Fourier pairs) to account for seasonality (see Section 2.4 of Kremser et al. (2014)), is fitted to the daily observational data to 5 Geosci. Model Dev. Discuss., doi:10.5194/gmd-2017-36, 2017 Manuscript under review for journal Geosci. Model Dev. Discussion started: 24 February 2017 c Author(s) 2017. CC-BY 3.0 License. obtain the mean annual cycle. The first 2 Fourier series expansions are given in Eq. 1.
where d is the day number of the year and β is a the regression coefficient being expanded. By using an offset basis function expanded in Fourier pairs, the resultant mean annual cycle is smooth. Examples of the mean annual cycle are shown in Figure   2 for four selected locations around New Zealand. This repeating mean annual cycle then provides the stationary baseline for the entire period of interest e.g. 1960 to 2100.

Training Phase
To determine the first-order long-term change in the magnitude of X, a simple linear correlation between X and T global is calculated for each of the 5 RCM simulations, viz.: which is used to calculate X , is generated using the same method and time period used to calculate the mean annual cycle of the observational set. X , rather than X, is used in Eq. 2 as the change in the seasonal cycle is of interest. Removing the mean annual cycle removed the need to add additional terms to describe the baseline seasonal cycle. The residuals from this fit are 10 then used by the stochastic weather generator (see Section 3.3) to model higher order changes in the variability in X which are not captured by Eq. (2).
Because GCM or RCM output provide a much longer time series than observations and extend into a period of greater changes in X, GCM or RCM output are preferentially used in this training phase. When using projections from RCMs as input to the training phase of EPIC, T global is sourced from the GCM that provided the boundary conditions for the RCM simulation.

15
Because the fit coefficient, α, is expected to depend on season, it is expanded in two Fourier pairs to account for its seasonality (Eq. 1). When embedded in Eq.
(2), this expansion results in five fit coefficients for α. We refer to each such set of five fit coefficients as a tuple; recall that this fit is applied at each grid point and for each available RCM simulation.
An example of a fit of Eq.
(2) to daily maximum surface temperature anomalies is shown in Figure 3 for a location in the South Island of New Zealand. The small annual cycle in the fit, with growing amplitude, results from summer-time and winter-time daily maximum surface temperatures exhibiting different correlations against T global . In addition to the long-term forced change, there is significant day-to-day variability. The use of the residuals from such fits in the stochastic weather generator model is described in Section 3.3.
The unitless α coefficient describes a locations sensitivity to changes in annual mean global mean surface temperature. The 5 magnitude of α indicates whether T max or T min are changing faster (α>1) or slower (α<1) than the global mean surface temperature. Example maps of the α coefficient for four selected days in New Zealand are shown in Figure 4. This analysis shows that daily maximum surface temperatures over most of New Zealand are warming slower than T global . However, high altitude regions, such as the Southern Alps, indicate T max increasing faster than T global for Spring, Summer and Autumn. There is, of course, some uncertainty in α. To account for that uncertainty, a large set of α tuples is derived through a Monte Carlo bootstrapping approach (Efron and Tibshirani , 1994), whereby residuals from the Eq. (2) fit are randomly sampled and added to the regression model fit to generate multiple statistically equivalent time series which are then refitted to obtain equally probable α fit coefficients (Bodeker and Kremser, 2015). This approach allows for the incorporation of the uncertainty in the fit of Eq. (2) into the final ensemble of projections.

5
It is also expected that α will depend on the RCM simulation used for the training. There are two ways in which this can be managed, viz.: 1. Tuples of α values are obtained for each RCM simulation providing data for the training phase of EPIC. Then, in the 'implementation phase' of EPIC (see below), for each ensemble member, a tuple of α values is selected from the Monte Carlo-derived set of α tuples from a RCM data set that is randomly selected.
2. Tuples of α values are obtained simultaneously for all RCMs by fitting Eq. (2) to concatenated time series of X and T global (y) obtained from all RCM simulations providing data for the training phase of EPIC.

5
For the purposes of this study, method (1) is used as method (2) will tend to underestimate the true uncertainty seen in the values of α.

Implementation Phase
Once the Monte Carlo derived sets (just one set if method (2) is used) of α tuples have been obtained, they are used in the implementation phase of EPIC. As described in Section 2.2, 190 simulations of T global can be generated using a SCM. A 10 randomly selected T global time series from the 190 member set is used together with a randomly selected tuple of α values to generate a series of maps of X f orced using Eq. (2), where the forced subscript denotes these are changes which correlate with T global There might be some concern that the random selection of an α tuple from the available set of tuples for a location could cause the spatial coherence in the forced signal across New Zealand to be lost, as at a nearby location a different tuple could be 15 randomly selected. This was tested for and was found not to be the case as the multiple instances of tuples (multiple instances of Figure 4) are very similar and consistent (not shown here).

Indirect response to T global and weather noise
In addition to the change in X that correlates directly with T global , higher order components of variability, as well as realistic weather noise, must be present in the projections comprising the ensemble. One potential use of the ensemble of projections 20 generated by EPIC is assessment of the impacts and implications of climate change on a regional scale. These impacts seldom happen at a single site, i.e. the impact is felt over a large area. For this reason it is important that any specific member of an ensemble is appropriately spatially coherent over multiple sites. This is not achieved if the method considers each site in isolation since any purely stochastically determined weather noise added to a site would not be spatially coherent at neighbouring sites. For this reason, an empirical orthogonal function (EOF) approach, described by Lorenz (1956), is used to reveal the 25 modes of variability inherent in a series of daily maps of X after dependence on T global has been removed, and their principal component (PC) time series. Hereafter we refer to these modes of variability as 'weather modes'. These weather modes, and stochastically generated PC time series (P C syn ), are then used to construct a stochastic weather generator which generates realistic weather noise. The following is recognised in the construction of the stochastic weather generator model: 1. That VCSN data will provide the most realistic representation of weather noise.
2. That RCM simulations will simulate how that weather noise is likely to evolve in response to climate change (represented by T global ).
3. That the RCM simulations will be imperfect in simulating the patterns of variability derived from the VCSN data. 4. That there will be weather modes whose amplitude and variability will respond to climate change as well as others which will not. We begin by conducting an EOF analysis on VCSN data that have been detrended by removing the first order trend and on residuals from the training phase in Section 3.2.1. Where the patterns of variability obtained from VCSN and RCM diverge is considered to be the cut-off point for where the RCM simulation can be taken to have any integrity with regard to simulating forced changes in weather noise. Visual inspection of the EOF maps derived from VCSN and RCM data suggested that the 10 first four modes of variability were well represented by the RCM simulations (see Figure 5).

11
Geosci. Model Dev. Discuss., doi:10.5194/gmd-2017Discuss., doi:10.5194/gmd- -36, 2017 Manuscript under review for journal Geosci. Model Dev. Discussion started: 24 February 2017 c Author(s) 2017. CC-BY 3.0 License. It is clear from Figure 5 that the RCM EOFs exhibit the same modes of weather variability as seen in the VCSN data up until EOF pattern 4. These first four patterns of variability together explain 83.3% of the total weather variability in the VCSN data and 64.7% of the variability in the RCM data. It is these four modes of weather variability that evolve with T global in our stochastic weather generator.

Modelling forced changes in the amplitude and variability of weather modes 5
To compare statistics from the PC time series calculated from VCSN and RCM data, they must share the same underlying weather modes. This is done by projecting the VCSN weather modes (EOF V CSN ) onto the RCM data to calculate a pseudo-PC time series. A pseudo-PC time series is calculated in the same way that a standard PC time series is calculated, except that the weather modes are prescribed instead of being calculated from the data. A pseudo-PC time series describes the magnitude of a particular pattern of variability from VCSN, which is present in the RCM data. The VCSN weather modes, rather than the 10 RCM weather modes, were prescribed because the observational data set is more likely representative of patterns of variability seen in New Zealand. The pseudo-PC and VCSN P C time series can be compared as they both describe the same patterns of variability.
In the stochastic weather generator, we consider changes in: 1. The amplitude of the weather mode: this is quantified by correlating the associated pseudo-PC time series with T global 15 and then using that correlation coefficient (β) to drive a trend in the PC time series obtained from the VCSN-based EOF analysis.
2. The variability of the weather mode: this is quantified by correlating the variability in the associated pseudo-PC time series with T global and then using that correlation coefficient (β var ) to drive a trend in the variability of the PC time series obtained from the VCSN-based EOF analysis. The mean variability of the weather mode is obtained from the VCSN PC 20 time series rather than the pseudo-PC time series, so that the weather mode emulates the magnitude of variability seen in the VCSN data.
We also recognize that the PC time series will exhibit temporal auto-correlation and therefore, that correlation is quantified and removed before correlating the PC signal, and its variability, against T global . The resulting time series (P C syn_n ) captures both long-term shifts and/or changes in spread of the n th weather mode. We note, however, that by considering only lag-one 25 autocorrelation in these PC time series, we neglect longer term auto-correlation, e.g. that resulting from El Niño and La Niña events. As a result, our ensemble time series exhibit smaller interannual variability than is observed in VCSN time series.
The ability of the method to generate a set of PDFs of the P C syn_1 to P C syn_4 time series is demonstrated in Figure 6. The EPIC method corrects for any shortcomings in the ability of the RCM to correctly simulate expected magnitudes of weather variability for these four primary modes and then accommodates these corrections when generating PC time series that evolve into the future.

Modelling higher order modes of variability in weather
The stochastic weather generator model includes the effects of EOF patterns five and higher but assumes that these modes 5 show no dependence on T global as the RCM simulations do not accurately simulate these higher modes of weather variability.
The variability of the PC time series often has a strong seasonal cycle. Therefore, for EOF pattern five and higher, synthetic PC time series (P C syn ) are generated using a standard Monte Carlo approach, i.e. randomly selecting values from N (0, σ(d)), that is, a normal distribution with a mean of 0 and a standard deviation which depends on the day of the year which is being modelled. σ(d) is determined by a zeroth-order fit expanded in two Fourier pairs to the VCSN PC time series. This approach allows selection of extreme PC values that are outside of the range of PC values experienced in the 1972-2013 period, but noting that the PDFs of these PCs do not evolve with time. As with the forced changes in the amplitude and variability of weather modes, the auto-correlation in the PC time series is also quantified and captured in the statistically modelled PC time 5 series.
For a given ensemble member, once synthetic PC time series at daily resolution have been generated, they are used to produce a reconstructed weather field, W , according to: where i, j and t represent the latitude, longitude and time dimensions respectively and n is the nth weather mode. 10 Since W has been constructed from a linear combination of spatial patterns of variability, each of which is spatially coherent, it retains the property of spatial coherence. The variability evolves as expected under changes in T global for the first four modes of variability, as simulated by the RCM, and where extreme conditions, outside the range of the training period, do occur with a statistically reasonable frequency due to the stochasticity in the construction of the pseudo PC time series.
T min is modelled identically to T max with one small change: days with anomalously low T max would be more likely to 15 have anomalously low T min . Not accounting for this correlation could result in stochastically modelled T min values being higher than the modelled T max value for that day. To avoid that, and to capture the correlation between T max and T min on any given day, the same set of random numbers used to generate the values in the synthetic PC n time series for T max for a given day is used to generate the values in the synthetic PC n time series for T min . This forces the selection of PC syn values from the same region of the PDF for both T max and T min . 20

Results
Examples of the T max and T min time series generated by the EPIC method are shown in Figure 7 for four population centres in New Zealand together with the associated VCSN time series. and therefore show no systematic bias with respect to the VCSN data. The EPIC-generated time series also show a long-term evolution consistent with expectations from RCM simulations, including the effects of the spread in those simulations. While it cannot be directly seen from the time series plotted in Figure 7, the EPIC-generated time series also exhibit changes in weather variability consistent with RCM projections of expected changes in the first four modes of weather variability. The apparent 5 annual cycle in the anomaly time series reflects the annual cycle in the variance and not an annual cycle in the anomalies; towards the end of the period there is a true annual cycle in the anomalies from differential seasonal changes in T max and T min . The interannual variability of the EPIC ensemble members is lower than that of the observational data set. This is due to EPIC not including any terms which describe patterns of variability which occur at time scales of longer than 1 year.

10
The EPIC (Ensemble Projections Incorporating Climate model uncertainty) method, is able to generate large ensembles of daily time series of daily maximum and minimum temperatures that exhibit the following characteristics: -No bias with respect to VCSN data.
-Long-term evolution consistent with projections from a suite of RCM simulations, incorporating the uncertainties inherent in those simulations as well as additional structural uncertainties that may arise from the use of a wider suite of 15 RCMs as captured by the use of projections of T global . T global time series were generated by a SCM tuned to 19 different AOGCMS and 10 different carbon cycle models and used as a predictor for the long-term change in T max and T min .
-Weather variability with extremes that extend beyond that observed in the VCSN record and that evolve in a way consistent with RCM projections of changes in the four primary modes of weather variability.
-Spatial coherence in weather variability in any single ensemble member is preserved. 20 As such, EPIC-generated projections are suitable for generating robust PDFs of projections of T max and T min .
The number of members in each ensemble is essentially limited only by the computing resources available. For calculating the PDFs that are delivered to users, we currently generate 19,000 members in each ensemble (one ensemble for each RCP scenario) at each 0.05 • × 0.05 • grid point across New Zealand.
A web-based tool has been developed to deliver PDFs of T max and T min for the period 2001-2010 and 2091-2100 to users 25 along with statistics regarding the change in frequency of extreme events, i.e. days per year with T max above 25 • C and 30 • C and T min below 0 • C and 2 • C. The tool is available at http://futureextremes.ccii.org.nz/.
The next steps for the development of EPIC include extending the range of climate variables to daily surface broadband radiation, surface humidity and precipitation, and incorporating longer term sources of variability, e.g. those generated by El Niño and La Niña events, into the stochastic weather model.