Fordham, D., Brown, S., Wigley, T., Rahbek, C., 2019. Cradles of diversity: unlikely relics of regional climate stability. Current Biology (https://doi.org/10.1016/j.cub.2019.04.001)
The R code necessary to re-produce the analysis is provided here.
Continuous paleoclimate simulations of gridded monthly-mean temperature at the surface (TS) for the period 21,000 BP to 100 BP were accessed from the TRaCE21ka experiment1 using PaleoView2. Multi-century pre-industrial control runs of annual-mean TS for 18 different AOGCMs (Table 1) were retrieved from the Coupled Model Intercomparison Project Phase 53 and resampled to a common 2.5° x 2.5° grid using bilinear interpolation. Bilinear interpolation was chosen because (i) the source and destination grids were rectilinear, (ii) temperature varies smoothly spatially, and (iii) bilinear interpolation reasonably retains the integrity and limitations of the original model output data, where orography is highly smoothed relative to the real-world4.
Pre-industrial control runs are multi-century unforced climate simulations, where the initial model conditions are set based on atmospheric gas concentrations pre-wide-scale industrialisation (c. 1750 - 1850 CE), combined with non-evolving boundary conditions3. We elected to use only the first realisation (r1i1p1) from each AOGCM for the pre-industrial control runs because all models, with the exception of the Community Climate System Model ver. 4 (CCSM4)5, only had a single realisation for pre-industrial control conditions. Furthermore, the additional pre-industrial realisations (r2i1p1 and r3i1p1) for the CCSM4 model ran for considerably shorter periods: 156 and 120 years (Table 1). These pre-industrial control runs were used as a baseline for determining centuries of large changes in global-mean temperature (see below).
We used multi-century pre-industrial control runs from the AOGCMs to estimate internally-generated natural variability in area-weighted global-mean trend magnitudes in TS(°C/year)6 using generalised least squares (GLS) regression and maximally (1-year) overlapping windows i.e., a 100-year moving window with a 1-year step between windows. GLS regression models were generated using the nlme package for R7. Trends for the 18 AOGCMs were calculated using area-weighted global-mean temperatures based on the slope of the GLS model. An AR(1) error structure was chosen to minimise any effect of temporal auto-correlation in the model residuals7. The resulting global pre-industrial model trends (i.e., the slope of the regression against time) for TS were used to generate a multi-model empirical cumulative distribution function (CDF) using unsigned (i.e. absolute) trends (Figure 1). Because the number of years varied between pre-industrial control runs from different AOGCMs (Table 1) we used a bootstrap procedure to ensure that all models had equal weights in the CDF. This involved randomly selecting 141 centuries from each model (with replacement) - 141 is the minimum number of maximally overlapping century-length windows across all models. The bootstrap procedure was repeated 1000 times. Centuries having extreme trends were defined as those where the trend exceeded the 90th percentile of the multi-model control-run CDF. The 90th percentile is commonly used in climate change research to identify extreme events8,9.
We used a similar approach to calculate paleoclimate trend magnitudes in annual TS (using the TraCE-21ka daily paleoclimate simulations) using maximally overlapping century "windows" spanning the period from 21,000 BP to 100 BP. The primary difference being that for the paleoclimate simulations we calculated 'local' measures of trend in annual mean surface temperature for each grid-cell (n = 10,368 cells) as well as the area-weighted global-mean trend magnitude (see Figure 2 for a schematic of the workflow). We also calculated grid-cell estimates of variability, where variability was defined as the standard deviation (SD) of the residuals about the local trend10. We did not bootstrap the paleoclimate data because we only had one model.
We applied the 90th percentile trend threshold of the multi-model control-run CDF to the time series of area-weighted global-mean paleoclimate trends (based on maximal overlapping windows) to identify specific centuries of extreme change in global-mean temperature during the period 21,000 BP to 100 BP. We then subset the time series of paleoclimate simulations based on these periods of extreme change in area-weighted global-mean temperature. In doing so, 'local' measures of trend and variability (i.e., at the grid cell) were constrained to centuries of unusually high trends in global-mean temperature.
To check that our definition of extreme centuries of change in global-mean temperature was appropriate, we determined the amount of time a millennium was considered 'extreme'. We did this by calculating the % of time that a millennium was characterised by trends ≥ 90th percentile of the control run trends (Figure 3). This confirmed that known large-scale climatic events during the last deglaciation (e.g. Bølling-Allerød, Younger Dryas, Holocene climate optimum etc.) were correctly identified as periods of extreme climate change in our analysis. We also examined the trend in global-mean temperature through time (Figure 3).
We calculated the median values of the century trends and variability about the trends for each grid cell for the periods of rapid change in global-mean temperature from 21,000 BP to 100 BP. We then mapped cells along a gradient of climate stability for land and ocean realms separately. The use of both trend and variability combines information about low frequency (trend) and high frequency (SD) variability. We classified cells according to eight levels of climate stability, ranging from 'extremely stable' to 'extremely unstable': (i) ≤ 10th percentile of both trend and variability; (ii) > 10th and ≤ 25th percentile of both trend and variability; (iii) > 25th and ≤ 50th percentile of both trend and variability; (iv) ≤ 50th percentile of trend and ≥ 50th percentile of variability; (v) ≥ 50th percentile of trend and ≤ 50th percentile of variability; (vi) ≥ 50th and < 75th percentile of both trend and variability; (vii) ≥ 75th and < 90th percentile of both trend and variability; (viii) ≥ 90th percentile of both trend and variability.
To test whether climate-related refugia are likely to have been stationary in time and space during the most recent deglaciation event, we calculated and mapped the percentage of time each grid cell was characterised by 'stable' and 'extremely stable' climatic conditions in the face of unusually high trends in global mean temperatures. If a window is centred on, say, year N, then we used the 4066 land and 6302 ocean gridded values of trend magnitude and SD to calculate the distributions of these data corresponding to this year, and determine the 10th and 25th percentile points for land and ocean biomes separately. For any given grid point, the climate was defined as 'stable' if both the trend magnitude and SD were ≤ the 25th percentile, and as 'extremely stable' if both were ≤ the 10th percentile. At any grid point then, for the 'stable' case, if 1 denotes stable and 0 denotes not stable, we have a time series of 0s and 1s that identify periods of stable climate. We calculated time spent in 'extremely stable' climate conditions in a similar way. We chose the 10th and 25th percentiles because they are commonly used in climate change research to identify extreme events9,11. However, there is no way of knowing with any certainty whether these thresholds relate directly to mechanisms responsible for spatiotemporal variation in biodiversity patterns. To determine how sensitive our results are to the choice of threshold we did a sensitivity analysis, whereby we varied the thresholds for stable and extremely stable climates ± 10% (n = 20 stratified samples), calculated the grid cell SD and mapped these values.
1 Liu, Z., Otto-Bliesner, B.L., He, F., Brady, E.C., Tomas, R., Clark, P.U., Carlson, A.E., Lynch-Stieglitz, J., Curry, W., Brook, E., et al. (2009). Transient Simulation of Last Deglaciation with a New Mechanism for Bølling-Allerød Warming. Science 325, 310-314.↩
2 Fordham, D.A., Saltré, F., Haythorne, S., Wigley, T.M.L., Otto-Bliesner, B.L., Chan, K., and Brook, B.W. (2017). PaleoView: A tool for generating continuous climate projections spanning the last 21,000 years at regional and global scales Ecography 40, 1348-1358.↩
3 Taylor, K.E., Stouffer, R.J., and Meehl, G.A. (2011). An Overview of CMIP5 and the Experiment Design. Bulletin of the American Meteorological Society 93, 485-498.↩
4 Fordham, D.A., Saltre, F., Haythorne, S., Wigley, T.M.L., Otto-Bliesner, B.L., Chan, K.C., and Brook, B.W. (2017). PaleoView: a tool for generating continuous climate projections spanning the last 21 000 years at regional and global scales. Ecography 40, 1348-1358.↩
5 Gent, P.R., Danabasoglu, G., Donner, L.J., Holland, M.M., Hunke, E.C., Jayne, S.R., Lawrence, D.M., Neale, R.B., Rasch, P.J., Vertenstein, M., et al. (2011). The Community Climate System Model Version 4. Journal of Climate 24, 4973-4991.↩
6 Gent, P.R., Danabasoglu, G., Donner, L.J., Holland, M.M., Hunke, E.C., Jayne, S.R., Lawrence, D.M., Neale, R.B., Rasch, P.J., Vertenstein, M., et al. (2011). The Community Climate System Model Version 4. Journal of Climate 24, 4973-4991.↩
7 Pinheiro J, Bates D, DebRoy S, Sarkar D, and R Core Team (2017). nlme: Linear and Nonliner Mixed Effects Models. R package version 3.1-131↩
8 Harris, R.M.B., Beaumont, L.J., Vance, T.R., Tozer, C.R., Remenyi, T.A., Perkins-Kirkpatrick, S.E., Mitchell, P.J., Nicotra, A.B., McGregor, S., Andrew, N.R., et al. (2018). Biological responses to the press and pulse of climate trends and extreme events. Nature Climate Change 8, 579-587↩
9 Zhang, X., Alexander, L., Hegerl, G.C., Jones, P., Tank, A.K., Peterson, T.C., Trewin, B., and Zwiers, F.W. (2011). Indices for monitoring changes in extremes based on daily temperature and precipitation data. WIREs Clim Change 2, 851-870↩
10 Nadeau, C.P., Urban, M.C., and Bridle, J.R. (2017). Coarse climate change projections for species living in a fine-scaled world. Global Change Biology 23, 12-24.↩
11 Hao, Z., AghaKouchak, A., and Phillips, T.J. (2013). Changes in concurrent monthly precipitation and temperature extremes. Environmental Research Letters 8, 034014↩