Abstract

The authors use a statistical regional climate model [Statistical Regional Model (STAR)] to project the Tibetan Plateau (TP) climate for the period 2015–50. Reanalysis datasets covering 1958–2001 are used as a substitute of observations and resampled by STAR to optimally fit prescribed linear temperature trends derived from the Max Planck Institute Earth System Model (MPI-ESM) simulations for phase 5 of the Coupled Model Intercomparison Project (CMIP5) under the representative concentration pathway 2.6 (RCP2.6) and RCP4.5 scenarios. To assess the related uncertainty, temperature trends from carefully selected best/worst ensemble members are considered. In addition, an extra projection is forced by observed temperature trends in 1958–2001. The following results are obtained: (i) Spatial average temperature will increase by 0.6°–0.9°C; the increase exceeds 1°C in all months except in boreal summer, thus indicating a reduced annual cycle; and daily minimum temperature rises faster than daily maximum temperature, resulting in a narrowing of the diurnal range of near-surface temperature. (ii) Precipitation increase mainly occurs in early summer and autumn possibly because of an earlier onset and later withdrawal of the Asian summer monsoon. (iii) Both frost and ice days decrease by 1–2 days in spring, early summer, and autumn, and the decrease of frost days on the annual course is inversely related to the precipitation increase. (iv) Degree-days increase all over the TP with peak amplitude in the Qaidam Basin and the southern TP periphery, which will result in distinct melting of the local seasonal frozen ground, and the annual temperature range will decrease with stronger amplitude in south TP.

1. Introduction

The Tibetan plateau (TP), with an average elevation of over 4000 m and an area of approximately 2.5 × 106 km2, exerts a huge influence on regional and global climate through thermal and mechanical forcing (Flohn 1968; Yeh and Gao 1979; Yanai et al. 1992; Ding 1992). It is underlain by plateau permafrost and has the largest cryospheric extent outside of the polar region supplying freshwater for all large Asian rivers, thus being recognized as the “third pole” and the “Asian water tower.” This setting marks the high sensitivity of this region to the ongoing warming (e.g., You et al. 2008a,b), which has been manifested through widespread permafrost warming and thawing and glacier melting in the last decades (Yao et al. 2007; Wu et al. 2013).

The response of these vulnerable elements of the climate system can further feed back to the climate system: for example, thawing permafrost will release the stored carbon into the atmosphere and thus enhance the ongoing warming (Anisimov 2007; Kennedy et al. 2008). Therefore, it is of great importance to assess the local expression of global warming in the TP and to obtain an insight into the possible warming extent in the near future, which should be taken into account for sustainable development planning.

So far this assessment is performed mainly based on numerical general circulation models (GCMs), which perform poorly in regions like the TP because of their coarse resolution, unresolved subgrid-scale processes, and the poor representation of the dramatic topographic gradients (Iorio et al. 2004). The poor performance of GCMs comes partly also from the delicacy of handling the glacier/snow cover–related albedo feedbacks in climate models (e.g., Qu and Hall 2007).

Another approach to understanding local climate variability is to nest a regional climate model (RCM) into a GCM; thus, a finer resolution can be achieved in the region of interest while the physical processes are also available for investigation (e.g., Ding et al. 2006; Gao et al. 2006; Ju et al. 2007; Srinivas et al. 2012). But although RCMs have shown improved performances in the simulation of present-day climate (Giorgi 2006; Laprise et al. 2008), noticeable biases have been noted even in their control simulations (e.g., Kotlarski et al. 2005) limiting their use for near-future regional climate impact studies, not to mention that small perturbations to boundary conditions can lead to significantly different future behavior (Collins and Allen 2002).

In comparison, statistical schemes are less exposed to these biases faced by GCMs and RCMs and emerge as a means complementary to numerical approaches. Examples of such statistical approaches are regression models (e.g., Stoner et al. 2013), weather pattern approaches (Wilby et al. 1998), and stochastic weather generators (Wilks 2010, 2012). The statistical resampling model used in this study [Statistical Regional Model (STAR)] resembles a weather generator in the sense that it resamples daily training data and generates ensembles of daily time series of temperature to optimally fit a linear trend prescribed for the projection period; other climate variables like precipitation from the training period are rearranged in the same order as the temperature time series, thus keeping physical consistency among all variables. It allows ensemble climate projections such that future climate statistics such as long-term averages or extreme events can be estimated. Since the resampling of past observations is conditioned on the prescribed temperature trend, STAR can be seen as an extended bootstrap approach (Efron and Tibshirani 1993). It has been successfully applied to the Elbe region and the Yangtze catchment (Orlowsky et al. 2008, 2010).

We use the dataset developed in the framework of the European Union Water and Global Change (WATCH) project. This WATCH dataset serves as a validation base for testing STAR's performance in the TP and also as a training archive for the future climate projections. Spatial features of the future climate are considered by applying STAR to multiple representative stations, and various climate variables are presented in terms of their time means, annual cycles, and spatial distributions. The article is organized as follows: section 2 briefly introduces the STAR model and the data used in this study; section 3 presents the validation results to demonstrate the applicability of STAR for the TP region; section 4 shows the observed present-day climate (training period 1958–2001) and the future projections (2015–50); section 5 focuses on discussing the results and caveats of the method; and section 6 presents the main conclusions.

2. Model and data

a. STAR

STAR was developed to generate regional climate projections for the near future by resampling the past climate time series. As shown in Fig. 1, given an observed time series of a characteristic climate variable, chosen as near-surface air temperature in our study, STAR assembles the observations in segments following two steps to optimally fit a linear trend prescribed for the future period: first, yearly segments (annual means, marked as X) are rearranged to obtain a first-order approximation of the linear trend (black line in the future period); second, the daily input time series are classified as blocks of certain length through a cluster analysis (12 days in our study) to keep synoptic information. STAR then works to improve the fitting by resampling these 12-day blocks (gray vertical lines) until one optimal fitting is obtained; other climate variables are subsequently assembled in the same order, thereby maintaining physical consistency among various climate variables. In this sense, generating a future series can be seen as defining a day-to-day mapping, via which the weather of any day in the future period is assigned the weather feature of a certain date from the observed series (Orlowsky et al. 2010).

Fig. 1.

STAR work chart at station 1 (location marked in red in Fig. 2a) zooming in on 1999–2001 and 2015–17. STAR first resamples (left) training data from 1958 to 2001 (black) in yearly segments (yearly means; marked as X) to fit (right) the prescribed temperature trend for 2015–50 (black line; with magnitude of 0.286°C decade−1) to a first order. Second, STAR resamples 12-day blocks (vertical light gray lines in both panels) to obtain a closer fitting to the linear trend. Step 2 is repeated until a given tolerance is fulfilled (in this case, within a distance of ±0.15°C from the given trend).

Fig. 1.

STAR work chart at station 1 (location marked in red in Fig. 2a) zooming in on 1999–2001 and 2015–17. STAR first resamples (left) training data from 1958 to 2001 (black) in yearly segments (yearly means; marked as X) to fit (right) the prescribed temperature trend for 2015–50 (black line; with magnitude of 0.286°C decade−1) to a first order. Second, STAR resamples 12-day blocks (vertical light gray lines in both panels) to obtain a closer fitting to the linear trend. Step 2 is repeated until a given tolerance is fulfilled (in this case, within a distance of ±0.15°C from the given trend).

In addition, a set of heuristic rules ensures that STAR simulated series, besides complying with the prescribed regression parameters (the linear temperature trends), exhibit realistic annual cycles and persistence. For example, a date window constrains the time span within which a selected date can be used as an analog such that weather on 1 January of a year can only possibly be assigned to a winter date (other than a summer date), as long as qualified dates are available, thus enforcing a realistic annual cycle. In the case of unavailability of such dates, weather blocks from other seasons are incrementally allowed, which often arises from overstretching of the training data.

The underlying assumption of STAR is that sequences of observations from the past can occur again in a very similar manner in the future: that is, the future climate and the observations share similar statistical features. This requirement is most likely not fulfilled in a drastically warmer climate, because samples (blocks) of warm weather will be more frequently used, thus leading to a biased sampling, which may lead to a narrower ensemble spread of temperature or even systematic cold biases in warm months in the case when available blocks are not warm enough. Yet, Orlowsky et al. (2010) point out that, as long as long-term statistics are considered, STAR yields satisfactory results if the warming in the training period continues to the near future with the same or even bigger strength. Their work serves as our guideline to determine which emission scenarios to pursue when applying STAR for future climate projections.

Theoretically, STAR can be applied to arbitrarily many stations, but the related computation cost rises tremendously with an increasing number of stations, since each optimal fitting has to satisfy all regression lines within reasonable tolerances (0.1°–0.15°C in our study), if it exits at all. In practice, STAR considers a limited number of stations to take into account spatial variability of the climate variables while keeping the computation cost in an acceptable range. These stations should, in essence, represent a sufficiently high portion of the total variability of the data in the region of interest. Various methods like the k-means cluster analysis method, implemented in the STAR package, or any other method that can extract centers of action of the underlying structure in the data can be used to determine these representative stations. Individual temperature trends are then prescribed to each representative station, and STAR searches for a resampling order that satisfies all prescribed trends. Since there are often more than one optimal fitting to the given regression line(s), STAR generates a spread of fitting by introducing stochastic elements thereby providing simulation ensembles. In this manner, STAR provides an efficient low-cost approach to climate uncertainty assessment. It is worth stressing that, besides determining the representative stations and their trends, no adaptation specific to the region of interest is necessary. More details about the STAR model and its merits and limitations have been discussed in Orlowsky et al. (2008, 2010).

b. Data

Ideally, STAR was planned for projecting future climate by making full use of observations. However, observational data with a reasonable spatial coverage and daily sampling frequencies are not available in the TP. Thus, reanalysis datasets are used in our study as a substitute of observational data, and the terminology “station” refers to grid cells.

1) Reanalysis data

The reanalysis data used was produced by the European Union WATCH project based on the 40-yr European Centre Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40). To obtain the WATCH data, the ERA-40 data were subject to sequential interpolation to ½° resolution, elevation correction, and monthly scale adjustments based on Climatic Research Unit (CRU; for corrected temperature, diurnal temperature range, and cloud cover) and Global Precipitation Climatology Centre (GPCC; for precipitation) monthly observations. Corrections for varying atmospheric aerosol loading and separate precipitation gauge for rainfall and snowfall have also been taken into account. Comparison with flux network (FLUXNET) data demonstrates a close correspondence between field-measured and these adjusted WATCH datasets for all variables (Weedon et al. 2010).

The WATCH data are used for both validation and future projections. In the validation experiment, WATCH data are divided into two periods, 1958–80 and 1981–2001. The former period is input to STAR, while the latter provides linear temperature trends and serves as observations to which STAR's performance is compared and thereby validated. For future projections, WATCH data from the entire period 1958–2001 serve as training data to fit linear temperature trends for 2015–50 derived from numerical models.

2) General circulation model data

We analyze simulations with the Max Planck Institute Earth System Model, low resolution (MPI-ESM-LR), devoted to phase 5 of the Coupled Model Intercomparison Project (CMIP5). MPI-ESM includes the atmospheric model ECHAM6 in T63 (1.9° × 1.9°) resolution with 47 vertical levels described by Stevens et al. (2013), the Max Planck Institute Ocean Model (MPI-OM) at approximately 1.6° resolution with 40 vertical layers (Jungclaus et al. 2006), and the land surface model Jena Scheme for Biosphere–Atmosphere Coupling in Hamburg (JSBACH; Raddatz et al. 2007) sharing the same horizontal grid as ECHAM6. All MPI-ESM modules interact directly without flux adjustments. Note that biogeochemical components are coupled, and thus the MPI-ESM simulations have the capability of using time-evolving emissions of constituents from which concentrations can be computed interactively. A detailed description of the model and an evaluation of the model performance regarding temperature and precipitation fields are given by Giorgetta et al. (2012). For an overview of CMIP5, readers are referred to Taylor et al. (2009, 2012).

Four groups of MPI-ESM simulations are analyzed: historical (HIST) and decadal prediction (DECAD) simulations and two scenarios, representative concentration pathway 2.6 (RCP2.6) and RCP4.5. A few details are given as follows: (i) HIST and DECAD simulated temperatures (spatial averages over the TP) for 1958–2001 are compared to the WATCH data to give a brief overview of their performances in the TP (Fig. 3). Note that MPI-ESM DECAD experiments are initialized every year with observed climate states for the period 1960–2001 and produce a set of 10-yr hindcasts (Fig. 3 presents only the 1-yr forecasts from the DECAD experiments). (ii) RCP2.6 and RCP4.5 provide temperature trends for 2015–50 that are prescribed to STAR for projection purposes. Details about how the temperature trends are estimated are introduced in section 2c.

3) Asian summer monsoon index

The daily Indian summer monsoon (ISM) data from 1958 to 2001 are downloaded from the Meteorological Department of the University of Hawaii at Manoa (Wang and Fan 1999; Wang et al. 2001) to describe the Asian summer monsoon. It is defined as the 850-hPa zonal wind shear between two regions: 5°–15°N, 40°–80°E and 20°–30°N, 70°–90°E.

c. Estimation of temperature trends

In our study, temperature trends are equivalent to temperature changes (difference between the end and start year of the linear regression line). Following the delta change approach (Hay et al. 2000), these temperature trends are added onto the mean temperature estimated from the last 15 yr of the WATCH data (1987–2001) and then prescribed to STAR for the projection period. By taking the 15-yr averages, we wish to reduce the influence of interannual variability. Each future projection has 100 ensembles.

To obtain realistic future climate projections, most relevant for STAR is a good estimation of future temperature trends. We take a practical approach: that is, to select for each representative station the ensemble member with the best performance in simulating linear temperature trends under the current climate, which we suppose will keep its good performance into the projection period and provide a reasonable projection of future temperature trends. To assess the related uncertainty, we select both the best and the worst ensemble members for each station, thereby bounding the unknown “true” future climate state: that is, the temperature evolution in this context.

Details of the trend estimation are given as follows: Temperature trends of three HIST ensembles are ranked stationwise according to their distance from the WATCH data in 1958–2001. The ensemble with the smallest/largest deviation is labeled as best/worst, forming two lists of five-member ensemble combinations shown as “best path” and “worst path” in Table 1. In the next step, we perform both best and worst projections with STAR for both RCP2.6 and RCP4.5 emission scenarios for the period 2015–50, with temperature trends taken from the ensembles on the two lists (Table 1). Thereby, we obtain four groups of projections: RCP2.6-BEST, RCP2.6-WORST, RCP4.5-BEST, and RCP4.5-WORST. In addition, we perform another experiment RCP-WATCH, in which observed temperature trends from the WATCH data (1958–2001) are prescribed to STAR to explore how the climate is like if the current warming rate continues with the same strength in the projection period.

Table 1.

Best and worst simulated ensemble member at each station, determined from their relative distance from the observed trend (WATCH data) for 1958–2001. This combination of ensemble members is kept for the future projection period (2015–50) for extracting future trends.

Best and worst simulated ensemble member at each station, determined from their relative distance from the observed trend (WATCH data) for 1958–2001. This combination of ensemble members is kept for the future projection period (2015–50) for extracting future trends.
Best and worst simulated ensemble member at each station, determined from their relative distance from the observed trend (WATCH data) for 1958–2001. This combination of ensemble members is kept for the future projection period (2015–50) for extracting future trends.

3. Model validation

The idea of the subsequent validation is to let STAR project the climate for the validation period 1981–2001 by fitting linear trends derived from the “reality”: the WATCH data from the same period. STAR reassembles climate segments of the independent preceding time span 1958–80 so that by comparing its projections to the WATCH data we can validate STAR's performance.

The following four variables are considered in this study: daily precipitation and daily maximum, mean, and minimum near-surface air temperature. We have extended our results by including derived climate variables, such as accumulated degree-days, defined as the yearly summation of mean temperature above zero, to characterize warming impacts on snow, glacier, and permafrost melting. Only grid points above 2500 m above mean sea level are analyzed.

In the following, we first introduce how multiple representative stations are selected and then briefly present validation results in terms of temperature and precipitation.

a. Representative stations

The TP is influenced by various dynamical climate systems: It is sandwiched between the subtropical westerly jet in the north and the tropical easterly jet in the south: both exhibit seasonal north–south shift in association with the reversal of strong Asian monsoon. In summer, abundant moisture supply is brought to the TP from the Arabian Sea and the Bay of Bengal by the Indian summer monsoon and from the South China Sea and the western Pacific by the East Asian summer monsoon. Besides, the midlatitude westerlies supply the northern part of the TP with moisture mostly originated from the North Atlantic (Zhu et al. 2010). In winter, the zonal orientation of the Himalayas blocks the synoptic-scale exchanges of warm tropical air with cold polar air centered near Siberia; the only avenue of air exchange is east of the Himalayas over Southeast Asia. These geographically varying influences from various dynamical systems lead to large spatial differences in regional climates in the TP, which we can take into account by introducing multiple stations into the STAR model. We apply the standard empirical orthogonal function (EOF) analysis to monthly surface temperature (1958–80) after removing the climatological annual cycle. EOF1 is characterized as a monopole with the center of action in central TP; both EOF2 and EOF3 show a dipole pattern oriented in the northwest–southeast and northeast–southwest directions, respectively. These three EOFs accounting for 82% of the total variance (i.e., 46%, 25%, and 11%, respectively) correspond to five centers of action, which have the highest loadings in the correlation map between the respective principal component (PC) and the temperature field and are then chosen as the representative stations (Figs. 2a–c).

Fig. 2.

Correlation maps between monthly surface temperature and its leading three PCs from the EOF analysis of the validation WATCH data (1981–2001): (a) PC-1, (b) PC-2, and (c) PC-3. Representative stations are chosen as correlation extremes and marked by red points; they are numbered as 1–5 and referred to in Table 1. The leading three EOFs account for 82% of the total variance (46%, 25%, and 11%, respectively); annual cycles are removed prior to the EOF analysis. The brown color background is the orography from MATLAB.

Fig. 2.

Correlation maps between monthly surface temperature and its leading three PCs from the EOF analysis of the validation WATCH data (1981–2001): (a) PC-1, (b) PC-2, and (c) PC-3. Representative stations are chosen as correlation extremes and marked by red points; they are numbered as 1–5 and referred to in Table 1. The leading three EOFs account for 82% of the total variance (46%, 25%, and 11%, respectively); annual cycles are removed prior to the EOF analysis. The brown color background is the orography from MATLAB.

The validation has been repeated with different sets of representative stations to test how sensitive STAR results are with respect to the number of stations and their locations, which may suffer from limitations of the EOF analysis (Richman 1986). The following results are obtained: (i) Considering more representative stations improves spatial details of the projections. However, the improvement beyond the first three leading EOFs (corresponding to five stations marked in Fig. 2) becomes marginal, whereas the related computation cost grows drastically. (ii) The influence of the exact locations is trivial. Five centroids determined by the k-means cluster analysis method or by the EOF method with the WATCH data covering 1958–2001 are similarly distributed as in Fig. 2: namely, one is located in central TP and the other four are located at the four corners of the plateau. Validation projections with these five centroids or with stations neighboring those marked in Fig. 2 or with those determined from yearly instead of monthly surface temperature data from the entire period (1958–2001) render results similar to those presented in Figs. 35.

Fig. 3.

Temperature and precipitation averaged over the TP (above 2500 m MSL) for the validation period 1981–2001. (a) Annual means (in total 21 yr) from the validation WATCH data (asterisks), STAR projections (gray), and MPI-ESM HIST (open circles) and DECAD (filled circles) simulations. Red, green, blue denote three ensembles of MPI-ESM simulations. (b) Deviations of STAR projections from the validation WATCH data in terms of 21-yr mean. Residual quantile–quantile plots for daily (c) temperature (°C) and (d) precipitation (mm day−1). The asterisk in (b) denotes STAR ensemble median.

Fig. 3.

Temperature and precipitation averaged over the TP (above 2500 m MSL) for the validation period 1981–2001. (a) Annual means (in total 21 yr) from the validation WATCH data (asterisks), STAR projections (gray), and MPI-ESM HIST (open circles) and DECAD (filled circles) simulations. Red, green, blue denote three ensembles of MPI-ESM simulations. (b) Deviations of STAR projections from the validation WATCH data in terms of 21-yr mean. Residual quantile–quantile plots for daily (c) temperature (°C) and (d) precipitation (mm day−1). The asterisk in (b) denotes STAR ensemble median.

Fig. 4.

(top) STAR projected (gray) and observed (red) annual cycles of temperature (°C) and precipitation (mm day−1) for 1981–2001. Dashed and black solid lines denote interquartile and ensemble medians, respectively. (bottom) STAR projections' deviations from the validation WATCH data.

Fig. 4.

(top) STAR projected (gray) and observed (red) annual cycles of temperature (°C) and precipitation (mm day−1) for 1981–2001. Dashed and black solid lines denote interquartile and ensemble medians, respectively. (bottom) STAR projections' deviations from the validation WATCH data.

Fig. 5.

Spatial distributions of (a),(c) temperature (°C) and (b),(d) precipitation (mm day−1) in the validation period 1981–2001. (a),(b) Climatological means from the validation WATCH data. (c),(d) Deviation of STAR ensemble median from the respective observed climatological mean in (a),(b). Shading in (c),(d) marks areas where STAR ensemble median represents observed values at the 5% significance level and are thus not statistically different. Thick gray lines mark 2500 m above mean sea level.

Fig. 5.

Spatial distributions of (a),(c) temperature (°C) and (b),(d) precipitation (mm day−1) in the validation period 1981–2001. (a),(b) Climatological means from the validation WATCH data. (c),(d) Deviation of STAR ensemble median from the respective observed climatological mean in (a),(b). Shading in (c),(d) marks areas where STAR ensemble median represents observed values at the 5% significance level and are thus not statistically different. Thick gray lines mark 2500 m above mean sea level.

b. Validation results

Being conditioned on the observed temperature trends from 1981 to 2001, STAR projected temperature and precipitation area averages well spread over the range of the validation WATCH data, whereas MPI-ESM simulations show systematic deviations from observations (Fig. 3a). STAR has a bias of around 0.1°C in temperature and less than 1% in precipitation (Fig. 3b). Figures 3c,d display residual quantile–quantile plots of spatially averaged temperature and precipitation, which are conventional quantile–quantile plots rotated clockwise by 45° around the diagonal line (Marzban et al. 2011). The quantile–quantile plots show large deviations for low temperature extremes and high precipitation events but otherwise remain flat enclosing zero (Figs. 3c,d), suggesting STAR's faithful representations of climatic means and its limited performance in climate extremes. The annual cycles of area average temperature and precipitation are well reproduced (Fig. 4). It is interesting to note that temperature is underestimated in HIST by 2°–3°C but overestimated in DECAD by 1°–2°C, while both experiments overestimate precipitation over the TP with DECAD performing slightly better (Fig. 3a).

Considering that only five representative stations are taken into account, STAR captures satisfactorily the northward decreasing feature of annual mean temperature and precipitation (Fig. 5), which are not statistically different from their observed values at the 5% significance level, thus sharing similar climatic means with the observations. It performs as good in projecting seasonal means of temperature but not that of precipitation (not shown). In summary, the results above reveal that, STAR can faithfully project climatic means in the TP, which gives us the confidence to apply it to projecting the mean climate for the near future.

4. Present and future climate

To derive temperature trends from MPI-ESM projection integrations, model grids falling inside a circle centered at each representative station within a radius of 1.25° are selected (1–2 grid cells) and linear trends are estimated from temperature averages across these selected grids. Trends from all three ensembles of both RCP2.6 and RCP4.5 scenario projections are shown in Figs. 6b,c. Trends are then taken from the best/worst ensembles following the two lists (Table 1) for each scenario projection. In addition, the scenario RCP-WATCH is forced by the observed trends in 1958–2001.

Fig. 6.

MPI-ESM simulated temperature trends at five representative stations: (a) historical simulations for 1958–2001 and (b) RCP2.6 and (c) RCP4.5 scenarios for 2015–50. Temperature trends of the training WATCH data (1958–2001) are shown in (a) as references; locations of representative stations are shown in Fig. 2.

Fig. 6.

MPI-ESM simulated temperature trends at five representative stations: (a) historical simulations for 1958–2001 and (b) RCP2.6 and (c) RCP4.5 scenarios for 2015–50. Temperature trends of the training WATCH data (1958–2001) are shown in (a) as references; locations of representative stations are shown in Fig. 2.

Main results are summarized as follows:

a. Spatial averages: Time means (Fig. 7)

1) Temperature and precipitation

All scenarios show that both temperature and precipitation will increase in the TP. The projected temperature increase (in terms of ensemble median) is within 0.57°–0.60°C in RCP2.6 scenarios and within 0.86°–0.93°C under RCP4.5 scenarios, while RCP-WATCH projects modest warming of around 0.78°C. It is interesting to note that, under each emission scenario, the temperature increase along the BEST path tends to be slightly lower than that along the WORST path, whereas the opposite holds for precipitation.

Fig. 7.

Changes of area averages from 2015 to 2050 from the reference WATCH data (1958–2001): (top) temperature and precipitation and (bottom) daily maximum and minimum temperature (Tmax and Tmin). Precipitation is shown in percentages. Asterisks denote STAR ensemble medians.

Fig. 7.

Changes of area averages from 2015 to 2050 from the reference WATCH data (1958–2001): (top) temperature and precipitation and (bottom) daily maximum and minimum temperature (Tmax and Tmin). Precipitation is shown in percentages. Asterisks denote STAR ensemble medians.

2) Daily minimum and maximum temperature (Tmin and Tmax)

Increase in Tmin is bigger than that in Tmax in all scenarios, thus indicating a reduced diurnal cycle. This difference is more evident in two WORST cases.

b. Spatial averages: Annual cycles (Fig. 8)

1) Temperature

Temperature increases by around 1°C through all months except in boreal summer when the increase is milder, suggesting a reduced annual cycle. We notice that the ensemble spread in July is considerably narrower than in other months, which relates to STAR's biased sampling, as will be discussed later.

Fig. 8.

Annual cycles of (a) temperature, (b) precipitation, (c) Indian summer monsoon index, (d) ice days, and (e) frost days for 2015–50. Panels (a),(b),(d),(e) denote the respective area averages. STAR projections are in gray (100 ensembles); dashed and black solid lines denote ensemble interquartile and medians, respectively. The training WATCH data (1958–2001) are shown in red as references. Bottom panels in (a)–(e) show deviations of STAR projections from the reference WATCH data.

Fig. 8.

Annual cycles of (a) temperature, (b) precipitation, (c) Indian summer monsoon index, (d) ice days, and (e) frost days for 2015–50. Panels (a),(b),(d),(e) denote the respective area averages. STAR projections are in gray (100 ensembles); dashed and black solid lines denote ensemble interquartile and medians, respectively. The training WATCH data (1958–2001) are shown in red as references. Bottom panels in (a)–(e) show deviations of STAR projections from the reference WATCH data.

2) Precipitation and ISM

Precipitation increases mainly in early summer and autumn, which coincides with a stronger ISM index. It suggests that the ISM may have an earlier onset and later withdrawal leading to the monsoon-related precipitation increase.

3) Ice days

Ice days are defined as days with daily maximum temperature below zero. Climatological ice days remain at zero from May through September and peak in January, exceeding 20 days. In the near future, ice days decrease in boreal winter months by 1–2 days with the strongest decrease in March and November.

4) Frost days

Frost days are defined as days when daily minimum temperature is lower than zero. Frost days occur almost every day in winter months and are lowest in summer months at around 10 days. The projected decrease by up to 2 days occurs in early summer and autumn, which is shared by all projections.

It is interesting to note that the decrease pattern of frost days on its annual cycle is almost an inversed picture of the increase pattern of precipitation. This inverse relationship has already been found in observations: Portmann et al. (2009) report a similar inverse relationship between trends in precipitation and temperature during May–June in the southern United States, which is more pronounced in daily maximum temperature and 90th percentile exceedance trends.

c. Spatial features: Time means (Fig. 9)

1) Temperature

All scenarios show spatially heterogeneous warming: both best scenarios show strongest warming in the northwest with amplitudes of 0.6°C in RCP2.6-BEST and 0.9°C in RCP4.5-BEST, whereas both worst scenarios show strongest warming in northeast TP, with amplitudes of 0.7°C in RCP2.6-WORST and 1°C in RCP4.5-WORST; RCP-WATCH shows strongest warming of 0.9°C in central TP.

Fig. 9.

Spatial distributions of temperature (°C), precipitation (%), ATR (°C), and degree-days. (top) Climatological means of the reference WATCH data for 1958–2001. Deviations of STAR projections from the respective WATCH climatological means for (top middle)–(bottom) RCP2.6-BEST, RCP2.6-WORST, RCP4.5-BEST, RCP4.5-WORST, and RCP-WATCH; precipitation deviations are in percentage with respect to its climatological mean in the top row. Shading in deviation plots denotes significant changes at the 5% level.

Fig. 9.

Spatial distributions of temperature (°C), precipitation (%), ATR (°C), and degree-days. (top) Climatological means of the reference WATCH data for 1958–2001. Deviations of STAR projections from the respective WATCH climatological means for (top middle)–(bottom) RCP2.6-BEST, RCP2.6-WORST, RCP4.5-BEST, RCP4.5-WORST, and RCP-WATCH; precipitation deviations are in percentage with respect to its climatological mean in the top row. Shading in deviation plots denotes significant changes at the 5% level.

2) Precipitation

All scenarios converge to the feature that rainfall has its strongest increase in central and east TP by 6%–10%. It is interesting to note that two best scenarios and RCP-WATCH suggest reduced precipitation in the western TP, which however does not occur in the two worst scenarios.

3) ATR

Annual temperature range (ATR) is the temperature difference between the warmest and coldest month. Climatological ATR is strongest in north TP and decreases southward. This pattern is opposite to that of the climatological precipitation field: that is, high ATR versus low precipitation. ATR decreases in all scenarios, suggesting a weakened annual cycle, consistent with the changes observed in the temperature annual cycle (Fig. 8). Its decrease has the strongest amplitude in the southeastern TP (around 0.4°–1.2°C).

4) Degree-days

A feature shared by all future projections is that degree-days will increase all over the TP with the strongest amplitude in the Qaidam Basin and the southern part of the TP.

5. Discussion

STAR's successful application on the Yangtze catchment (Orlowsky et al. 2010) has been extended to the crucial climate region of the Tibetan Plateau, which renders the following results:

a. Time means of spatial averages

The area-averaged temperature will increase by 0.6°–0.9°C; daily minimum temperature rises faster than daily maximum temperature; thus, the temperature diurnal cycle will weaken, which has already been found in observations from 1996 to 2002 (Oku et al. 2006) and, according to STAR's projections, will continue into the projection period.

b. Annual cycles of spatial averages

The temperature increases through the year with weaker amplitude in boreal summer, suggesting a reduced annual cycle; precipitation increase occurs mainly in early summer and autumn, likely resulting from earlier onset and later withdrawal of the Asian summer monsoon, as has been observed in recent decades (Sabeerali et al. 2011; Kajikawa et al. 2012) and is foreseen by multiple GCMs in a warmer climate (Kitoh 2006). On the other hand, it is important to note that the mean state of the Asian summer monsoon and its variability and possible change under future warming varies considerably depending on the model examined and the monsoon index analyzed (e.g., Bao 2012) and remains one of the challenging tasks for future investigations.

STAR projects a decrease in frost days mainly in early summer and autumn but a concurrent increase in precipitation. This inversed relationship has also been found in the last few decades (Zhai and Pan 2003; Duan and Wu 2006) and is likely due to the increase in daily minimum temperature (and thus reduced frost days), which results in increased low clouds at night that favors precipitation.

c. Spatial features of time means

All scenarios project increased precipitation in central and east TP, whereas in west TP precipitation decreases in RCP-WATCH and the two best scenarios but not in the two worst scenarios. Another convergent feature among different scenarios is a projected overall warming by 0.5°–1.2°C and a decrease in ATR by 0.4°–1.2°C. The strongest ATR decrease occurs in southeast TP where relatively lower climatological ATR values are found. This reduction is mainly due to stronger warming in cold seasons and milder warming in warm months: for example, in RCP4.5-BEST, temperature increases by 0.9°–1°C in December–February and by 0.4°–0.5°C in June–August (not shown). Similar changes are also projected by MPI-ESM simulations (not shown). Degree-days will increase over the entire TP with enhanced amplitude in the Qaidam Basin and south TP, where seasonal frozen ground is distributed (Cheng and Wu 2007). This suggests that the seasonal frozen ground area is likely to undergo distinct warming and thus melting in the forthcoming decades.

Four of the five scenarios performed in this study are conditioned by the performance of MPI-ESM simulations in the historical period and the near future. The rationale for taking this single-model approach, instead of a multimodel ensemble method, is that we wish to explore and demonstrate the potential of combining merits of both numerical and statistical climate models, rather than to characterize the full uncertainties in climate projections, for which the latter method may serve better.

Uncertainty in climate projections were taken into account in two ways: (i) via selecting the best/worst ensembles and thereby defining the boundary of the related uncertainty spread under a given emission scenario and (ii) via generating 100 ensembles for each scenario. These two ways of considering uncertainty can be easily practiced in inexpensive personal computers, which is a merit reserved for statistical downscaling models.

There are a few points to be clarified. STAR does not know specific characteristics of the region studied, such as topography. We would like to mention two points on this: first, the observation surrogate, the WATCH data, contains a smoothed version of real topography, which already alleviated the effort of dealing with topography-related complications. Second, we confined our analysis only to grid cells above 2500 m to reduce the topography impact. Actually, when all grid cells above 500 m are considered, large biases occur in the southern TP periphery (not shown). Therefore, we like to stress that the choice of 2500 m is a compromise between the spatial coverage we wish to resolve and the complication introduced by complex topography.

It is foreseeable that when high spatial resolution data are available in regions like the TP, topography-related complications and other subtle feedback processes, such as snow/glacier cover–related processes, will become increasingly important, which certainly poses a great challenge for simple statistical models (e.g., STAR) and equally so for complex numerical models. In fact, this concern raises the importance of a rigid model validation, which enables users to know with certain confidence, where the strength and weakness of the model lie, when it is applicable at all.

In the end, it is important to note a few limitations of STAR. The hypothesis behind STAR and other statistical downscaling methods is that the future climate will not be fundamentally different from that in the observations so that the statistical and physical relationships between climate variables established for the present climate will persist into the near future. Therefore, if the future climate does not share the statistical features of the training climate, STAR will not be suitable any longer for climate projection purposes and leads to overstretching of the training data, which limits the number of suitable blocks leading to a biased sampling. In a warmer-future scenario, this overstretching of the training data tends to occur first in warm months, as shown by the narrowing of ensemble spread in July (Fig. 8). In fact, the milder warming in July may be a manifestation of cold biases as a result of this caveat of STAR, which may partly contribute to the reduction of the temperature annual cycle. For the very same concern, we decided not to consider the extreme scenario RCP8.5 but the peak-and-decay and intermediate scenarios (RCP2.6 and RCP4.5) instead. The latter two scenarios project milder temperature increase over the TP until 2050, with amplitude comparable to that of the observed trends (Fig. 3), whereas RCP8.5 projects a temperature increase that doubles the observed trend from 1958 to 2001 (not shown) and thus has a climate statistically different from the observations.

Another evident limitation of STAR is that the current version considers only linear trends, whereas in reality Earth climates exert natural variability on decadal and longer time scales (e.g., Shen et al. 2011). This constraint in usage of STAR should be kept in mind, which confined us to consider only long-term mean climate statistics. On the other hand, it also breeds the possibility of implementing a more realistic trend, which is attractive to consider as a future project.

6. Conclusions

We have demonstrated how STAR makes regional climate projections by combining observational data and carefully estimated future temperature trends, which both can be specifically tailored for the region of interest. The following results are obtained for the period 2015–50: (i) Spatial average temperature will increase by 0.6°–0.9°C. This warming occurs through all months with its weakest amplitude in boreal summer, thus indicating a reduced annual cycle. Both daily minimum and maximum temperature will rise leading to a reduction in ice and frost days (by 1–2 days), respectively. Further, daily minimum temperature rises faster than daily maximum temperature; thus, the diurnal range of near-surface temperature will become narrower. (ii) On the annual course, the decrease in frost days is inversely related to precipitation increase. This inverse relationship may be due to the increase in daily minimum temperature that favors low clouds at night and consequently more precipitation. On the other hand, the increased precipitation in early summer and autumn is primarily due to an earlier onset and later withdrawal of the Asian summer monsoon. (iii) Degree-days increase all over the TP with peak amplitude in the Qaidam Basin and the southern TP periphery, which may result in distinct melting of the local seasonal frozen ground. The annual temperature range will decrease with stronger amplitude in the south TP.

Acknowledgments

Wang is supported by the “Hundred Talent Program” of Chinese Academy of Sciences. Thanks to Hermann Österle for providing the data analyzed and to Frank Sielmann for helpful discussions. The authors would like to express thanks to anonymous reviewers for their helpful suggestions and comments.

REFERENCES

REFERENCES
Anisimov
,
O. A.
,
2007
:
Potential feedback of thawing permafrost to the global climate system through methane emission
.
Environ. Res. Lett.
,
2
, 045016,
doi:10.1088/1748-9326/2/4/045016
.
Bao
,
Q.
,
2012
:
Projected changes in Asian summer monsoons in RCP scenarios of CMIP5
.
Atmos. Oceanic Sci. Lett.
,
5
,
43
48
.
Cheng
,
G.
, and
T.
Wu
,
2007
:
Response of permafrost to climate change and their environmental significance, Qinghai-Tibet Plateau
.
J. Geophys. Res.
,
112
, F02S03,
doi:10.1029/2006JF000631
.
Collins
,
M.
, and
M. R.
Allen
,
2002
:
Assessing the relative roles of initial and boundary conditions in interannual to decadal climate predictability
.
J. Climate
,
15
,
3104
3109
.
Ding
,
Y.
,
1992
:
Effects of Qinghai-Xizang (Tibetan) Plateau on the circulation features over the plateau and its surrounding areas
.
Adv. Atmos. Sci.
,
9
,
112
130
.
Ding
,
Y.
, and
Coauthors
,
2006
:
Multi-year simulations and experimental seasonal predictions for rainy seasons in China by using a nested regional climate model (RegCM-NCC). Part I: Sensitivity study
.
Adv. Atmos. Sci.
,
23
,
323
341
.
Duan
,
A.
, and
G.
Wu
,
2006
:
Change of cloud amount and the climate warming on the Tibetan Plateau
.
Geophys. Res. Lett.
,
33
,
L22704
,
doi:10.1029/2006GL027946
.
Efron
,
B.
, and
R. J.
Tibshirani
,
1993
:
An Introduction to the Bootstrap. Monogr. Stat. Appl. Probab.,
No. 57, Chapman & Hall, 436 pp.
Flohn
,
H.
,
1968
:
Contributions to a meteorology of the Tibetan highlands
. Colorado State University Atmospheric Science Paper 130, 120 pp.
Gao
,
X.
,
J. S.
Pal
, and
F.
Giorgi
,
2006
:
Projected changes in mean and extreme precipitation over the Mediterranean region from a high resolution double nested RCM simulation
.
Geophys. Res. Lett.
,
33
, L03706,
doi:10.1029/2005GL024954
.
Giorgetta
,
M. A.
, and
Coauthors
,
2013
: Climate and carbon cycle changes from 1850 to 2100 in MPI-ESM simulations for the Coupled Model Intercomparison Project phase 5. J. Adv. Model. Earth Syst., doi:10.1002/jame.20038, in press.
Giorgi
,
F.
,
2006
:
Regional climate modeling: Status and perspectives
.
J. Phys.
,
139
,
101
118
.
Hay
,
L. E.
,
R. L.
Wilby
, and
G. H.
Leavesley
,
2000
:
A comparison of deltal change and downscaled GCM scenarios for three mountainous basins in the United States
.
J. Amer. Water Resour. Assoc.
,
36
,
387
397
.
Iorio
,
J. P.
,
P. B.
Duffy
,
B.
Govindasamy
,
S. L.
Thompson
,
M.
Khairoutdinov
, and
D.
Randall
,
2004
:
Effects of model resolution and subgrid-scale physics on the simulation of precipitation in the continental United States
.
Climate Dyn.
,
23
,
243
258
.
Ju
,
L. X.
,
H. K.
Wang
, and
D. B.
Jiang
,
2007
:
Simulation of the Last Glacial Maximum climate over East Asia with a regional climate model nested in a general circulation model
.
Palaeogeogr. Palaeoclimatol. Palaeoecol.
,
248
,
376
390
.
Jungclaus
,
J. H.
, and
Coauthors
,
2006
:
Ocean circulation and tropical variability in the coupled model ECHAM5/MPI-OM
.
J. Climate
,
19
,
3952
3972
.
Kajikawa
,
Y.
,
T.
Yasunari
,
S.
Yoshida
, and
H.
Fujinami
,
2012
:
Advanced Asian summer monsoon onset in recent decades
.
Geophys. Res. Lett.
,
39
,
L03803
,
doi:10.1029/2011GL050540
.
Kennedy
,
M.
,
D.
Mrofka
, and
C.
von der Borch
,
2008
:
Snowball Earth termination by destabilization of equatorial permafrost methane clathrate
.
Nature
,
453
,
642
645
.
Kitoh
,
A.
,
2006
: Asian monsoons in future. The Asian Monsoon, B. Wang, Ed., Praxis, 631–649.
Kotlarski
,
S.
,
A.
Block
,
U.
Böhm
,
D.
Jacob
,
K.
Keuler
,
R.
Knoche
,
D.
Rechid
, and
A.
Walter
,
2005
:
Regional climate model simulations as input for hydrological applications: Evaluation of uncertainties
.
Adv. Geosci.
,
5
,
119
125
.
Laprise
,
R.
, and
Coauthors
,
2008
:
Challenging some tenets of regional climate modelling
.
Meteor. Atmos. Phys.
,
100
,
3
22
,
doi:10.1007/s00703-008-0292-9
.
Marzban
,
C.
,
R.
Wang
,
F.
Kong
, and
S.
Leyton
,
2011
:
On the effect of correlations on rank histograms: Reliability of temperature and wind speed forecasts from finescale ensemble reforecasts
.
Mon. Wea. Rev.
,
139
,
295
310
.
Oku
,
Y.
,
H.
Ishikawa
,
S.
Haginoya
, and
Y.
Ma
,
2006
:
Recent trends in land surface temperature on the Tibetan Plateau
.
J. Climate
,
19
,
2995
3003
.
Orlowsky
,
B.
,
F.-W.
Gerstengarbe
, and
P.
Werner
,
2008
:
A resampling scheme for regional climate simulations and its performance compared to a dynamical RCM
.
Theor. Appl. Climatol.
,
92
,
209
223
.
Orlowsky
,
B.
,
O.
Bothe
,
K.
Fraedrich
,
F.-W.
Gerstengarbe
, and
X.
Zhu
,
2010
:
Future climates from bias-bootstrapped weather analogues: An application to the Yangtze River basin
.
J. Climate
,
23
,
3509
3524
.
Portmann
,
R. W.
,
S.
Solomon
, and
G. C.
Hegel
,
2009
:
Spatial and seasonal patterns in climate change, temperatures, and precipitation across the United States
.
Proc. Natl. Acad. Sci. USA
,
106
,
7324
7329
.
Qu
,
X.
, and
A.
Hall
,
2007
:
What controls the strength of snow-albedo feedback?
J. Climate
,
20
,
3971
3981
.
Raddatz
,
T. J.
, and
Coauthors
,
2007
:
Will the tropical land biosphere dominate the climate–carbon cycle feedback during the twenty-first century?
Climate Dyn.
,
29
,
565
574
.
Richman
,
M. B.
,
1986
:
Rotation of principal components
.
J. Climatol.
,
6
,
293
335
.
Sabeerali
,
C. T.
,
S. A.
Rao
,
R. S.
Ajayamohan
, and
R.
Murtugudde
,
2011
:
On the relationship between Indian summer monsoon withdrawal and Indo-Pacific SST anomalies before and after 1976/1977 climate shift
.
Climate Dyn.
,
39
,
841
859
,
doi:10.1007/s00382-011-1269-9
.
Shen
,
C.
,
W.-C.
Wang
, and
G.
Zeng
,
2011
:
Decadal variability in snow cover over the Tibetan Plateau during the last two centuries
.
Geophys. Res. Lett.
,
38
,
L10703
,
doi:10.1029/2011GL047288
.
Srinivas
,
C. V.
,
D.
Hariprasad
,
D. V.
Bhaskar Rao
,
Y.
Anjaneyulu
,
R.
Baskaran
, and
B.
Venkatraman
,
2012
:
Simulation of the Indian summer monsoon regional climate using advanced research WRF model
.
Int. J. Climatol.
,
33
,
1195
1210
,
doi:10.1002/joc.3505
.
Stevens
,
B.
, and
Coauthors
,
2013
:
Atmospheric component of the MPI-M Earth System Model: ECHAM6
.
J. Adv. Model. Earth Syst.
,
5
,
146
172
,
doi:10.1002/jame.20015
.
Stoner
,
A. M. K.
,
K.
Hayhoe
,
X.
Yang
, and
D. J.
Wuebbles
,
2013
:
An asynchronous regional regression model for statistical downscaling of daily climate variables
.
Int. J. Climatol.
,
33 (11), 2473–2494, doi:10.1002/joc.3603
.
Taylor
,
K. E.
,
R. J.
Stouffer
, and
G. A.
Meehl
,
2009
: A summary of the CMIP5 experiment design. PCDMI Rep., 33 pp.
Taylor
,
K. E.
,
R. J.
Stouffer
, and
G. A.
Meehl
,
2012
:
An overview of CMIP5 and the experiment design
.
Bull. Amer. Meteor. Soc.
,
93
,
485
498
.
Wang
,
B.
, and
Z.
Fan
,
1999
:
Choice of South Asian summer monsoon indices
.
Bull. Amer. Meteor. Soc.
,
80
,
629
638
.
Wang
,
B.
,
R.
Wu
, and
K. M.
Lau
,
2001
:
Interannual variability of the Asian summer monsoon: Contrast between the Indian and the western North Pacific–East Asian monsoons
.
J. Climate
,
14
,
4073
4090
.
Weedon
,
G. P.
,
S.
Gomes
,
P.
Viterbo
,
H.
Oesterle
,
J. C.
Adam
,
N.
Bellouin
,
O.
Boucher
, and
M.
Best
,
2010
: The WATCH forcing data 1958-2001: A meteorological forcing dataset for land surface- and hydrological-models. WATCH Tech. Rep. 22, 41 pp.
Wilby
,
R. L.
,
T. M. L.
Wigley
,
D.
Conway
,
P. D.
Jones
,
B. C.
Hewitson
,
J.
Main
, and
D. S.
Wilks
,
1998
:
Statistical downscaling of general circulation model output: A comparison of methods
.
Water Resour. Res.
,
34
,
2995
3008
.
Wilks
,
D. S.
,
2010
:
Use of stochastic weather generators for precipitation downscaling
.
Wiley Interdisc. Rev. Climate Change
,
1
,
898
907
,
doi:10.1002/wcc.85
.
Wilks
,
D. S.
,
2012
:
Stochastic weather generators for climate-change downscaling. Part II: Multivariable and spatially coherent multisite downscaling
.
Wiley Interdisc. Rev. Climate Change
,
3
,
267
278
,
doi:10.1002/wcc.167
.
Wu
,
T.
,
L.
Zhao
,
R.
Li
,
Q.
Wang
,
C.
Xie
, and
Q.
Pang
,
2013
:
Recent ground surface warming and its effects on permafrost on the central Qinghai-Tibet Plateau
.
Int. J. Climatol.
,
33
,
920
930
.
Yanai
,
M. H.
,
C.
Li
, and
Z.
Song
,
1992
:
Seasonal heating of the Tibetan Plateau and its effects on the evolution of the Asian summer monsoon
.
J. Meteor. Soc. Japan
,
70
,
189
221
.
Yao
,
T.
,
J.
Pu
,
A.
Lu
,
Y.
Wang
, and
W.
Yu
,
2007
:
Recent glacial retreat and its impact on hydrological processes on the Tibetan Plateau, China, and surrounding regions
.
Arct. Antarct. Alp. Res.
,
39
,
642
650
.
Yeh
,
T. C.
, and
Y. X.
Gao
,
1979
: Meteorology of the Qinghai-Xizang (Tibet) Plateau (in Chinese). Science Press, 278 pp.
You
,
Q. L.
,
S.
Kang
,
E.
Aguilar
, and
Y.
Yan
,
2008a
:
Changes in daily climate extremes in the eastern and central Tibetan Plateau during 1961–2005
.
J. Geophys. Res.
,
113
, D07101,
doi:10.1029/2007JD009389
.
You
,
Q. L.
,
S.
Kang
,
N.
Pepin
, and
Y.
Yan
,
2008b
:
Relationship between trends in temperature extremes and elevation in the eastern and central Tibetan Plateau, 1961–2005
.
Geophys. Res. Lett.
,
35
, L04704,
doi:10.1029/2007GL032669
.
Zhai
,
P.
, and
X.
Pan
,
2003
:
Trends in temperature extremes during 1951–1999 in China
.
Geophys. Res. Lett.
,
30
,
1913
,
doi:10.1029/2003GL018004
.
Zhu
,
X.
,
O.
Bothe
, and
K.
Fraedrich
,
2010
:
Summer atmospheric bridging between Europe and East Asia: Influences on drought and wetness on the Tibetan Plateau
.
Quat. Int.
,
236
,
151
157
,
doi:10.1016/j.quaint.2010.06.015
.