Multivariate climate departures have outpaced univariate changes across global lands

Changes in individual climate variables have been widely documented over the past century. However, assessments that consider changes in the collective interaction amongst multiple climate variables are relevant for understanding climate impacts on ecological and human systems yet are less well documented than univariate changes. We calculate annual multivariate climate departures during 1958–2017 relative to a baseline 1958–1987 period that account for covariance among four variables important to Earth’s biota and associated systems: annual climatic water deficit, annual evapotranspiration, average minimum temperature of the coldest month, and average maximum temperature of the warmest month. Results show positive trends in multivariate climate departures that were nearly three times that of univariate climate departures across global lands. Annual multivariate climate departures exceeded two standard deviations over the past decade for approximately 30% of global lands. Positive trends in climate departures over the last six decades were found to be primarily the result of changes in mean climate conditions consistent with the modeled effects of anthropogenic climate change rather than changes in variability. These results highlight the increasing novelty of annual climatic conditions viewed through a multivariate lens and suggest that changes in multivariate climate departures have generally outpaced univariate departures in recent decades.

Here, we examine trends in multivariate annual climate departures for global land surfaces during 1958-2017 relative to 1958-1987 baseline conditions. We focus on four variables that encapsulate the climatic basis for Earth's biota and have distinct impacts on human and natural systems: the coldest monthly average minimum temperature (T n,min ), warmest monthly average maximum temperature (T x,max ), annual actual evapotranspiration (AET), and annual climatic water deficit (D). These four variables collectively span thermal and moisture climatic axes and define the climatic niche of many species as well as climate impacts thereof. For instance, T n,min influences the life history and distributions of species due to the physiological, ecological, and evolutionary impacts 26 and is the basis of hardiness zones in the agricultural sector, while T x,max has been shown to have impacts on electricity demand 27 and crop yields 28 . AET and D are a reduced set of biologically relevant climate variables that represent the supply and unmet atmospheric demand components of the water balance 29 . AET is a proxy for productivity in natural systems 30 , while D is widely-used metric for ecosystem disturbance and the terrestrial carbon cycle 19,31 . In addition, various combinations of these variables are predictor variables for modeling species distributions 32,33 and climate impacts 34,35 .
We use Mahalanobis distance and its standardization based on the Chi distribution 3,36 to compute annual multivariate climate departures (σ d ). The Mahalanobis distance compactly quantifies multiple variables, accounts for their covariance structure, and measures the distance in multivariate space away from a centroid through principal components analysis of standardized anomalies. σ d can be interpreted as a multivariate z-score which scales climate departure with respect to local interannual climate variability. We compare σ d to its univariate counterpart for individual variables using standardized Euclidean distance (e.g., σ Tx,max ). To quantify the relative contribution of anthropogenic climate forcing to trends in climate departures, we develop a counterfactual simulation that removes the modeled climate change signal from observations. Two sensitivity experiments are used to elucidate the relative contribution of changes in means and changes in intra-annual to inter-annual variability on observed changes in climate departures. Both experiments remove the linear trends for variables during 1958-2017; the first experiment excludes linear trends in annual mean temperature, precipitation, and reference evapotranspiration, while the second only excludes linear trends in annual precipitation. Lastly, we compare observed trends in climate departures with those from internal climate variability using a 500 year control climate simulation.

Results
Positive trends in climate departures are evident in recent decades, particularly for σ d (Fig. 1a,b). The four individual calendar years with the greatest global median σ d coincided with the four warmest years globally over the period of analysis (2010,(2015)(2016)(2017). The 60-year median trend in σ d across global terrestrial surfaces was +0.78σ. Positive trends in climate departures over the 60-year period were observed for T x,max , T n,min , AET and D, but the magnitude of global median trends in σ d was approximately three times greater than trends in climate departures for individual variables (Fig. 1e). Similarly, the geographic extent of land with σ d exceeding 2 standard deviations increased markedly from a baseline of ~5% of land during the reference period  to ~30% during the most recent decade (2008-2017); increases in the geographic coverage of land >2 standard deviations for σ d outpaced increases for individual variables (Fig. 1c,d,f). Similar results were seen using a truncated non-parametric approach for transforming individual variables, although the magnitude of trends was reduced due to the conservative approach (Supplemental Methods). In addition to observed trends, we found a weak correlation (r = 0.21, p = 0.1) between global median σ d and the Multivariate ENSO Index 37 averaged over Jan-Jun, highlighting the tendency for larger climate departures during El Niño years.
Widespread positive trends in σ d were observed over the past six decades across most terrestrial surfaces with statistically significant increases present for 58% of lands and the largest increases in southern Europe, southeast Asia, Africa, and the Amazon ( Fig. 2a; Supplementary Fig. 1). Significant positive trends in σ Tx,max and σ Tn,min were seen across 16-18% of global lands (Fig. 2b,c). Significant positive trends in σ AET were seen for 16% of global lands, primarily across high latitudes of the Northern Hemisphere and equatorial regions. Significant positive trends in σ D were found for 11% of global lands. In contrast, significant negative trends in climate departures were confined to small geographic areas of 1-5% of land area depending on variable.
Geographic hotspots of large positive trends often occurred in areas with low interannual variability. This is further demonstrated by examining trends in climate departures and the magnitude of interannual variability by latitude and biome. Latitudinal patterns show large positive trends in σ d and σ Tn,min near the equator co-located with regions with low variability during the reference period (Fig. 3a). Additionally, larger positive σ d trends were seen in the Northern Hemisphere than the Southern Hemisphere (Fig. 3a). Trends for temperature-based departures generally exceeded moisture-based departures from 40°S to 50°N. At higher-latitudes, trends in σ AET outpaced trends in temperature departures consistent with higher temperature variability at higher latitudes (Fig. 3b). Trends in climate departures by biome were generally largest where interannual variability among the four climate variables was low (Fig. 3c,d), with the largest positive σ d trend in the Tropical Forest biome and the smallest positive trend in the Temperate Grassland biome.
Larger positive σ d trends compared to univariate climate departure trends can be understood by considering the covariance structure and trends in the underlying climate variables. AET and D exhibited strong negative correlations on interannual timescales across most terrestrial surfaces, with T x,max exhibiting positive correlation to D ( Supplementary Fig. 3). By contrast, T n,min was weakly correlated to the other climate variables, consistent with T n,min being seasonally out of phase with the primary climatic factors that influence surface water availability. The principal components (PC) used in the Mahalanobis distance calculation largely reflect these relationships with the first loading projecting onto patterns of water limitations, loadings 2 and 3 projecting onto orthogonal dimensions related to T n,min and T x,max , respectively, and loading 4 projecting to same signed anomalies of AET and D ( Supplementary Fig. 4). Positive trends in T x,max and T n,min were seen across most terrestrial surfaces with AET and D trends showing more spatial variability but generally positive trends during 1958-2017 ( Supplementary  Fig. 5). Trends in σ d , by definition, arise from trends in Euclidean distances of PC scores. While no single loading predominantly accounted for σ d trends, trends in PC4 and PC3 accounted for a majority of positive σ d trends for 22% and 16% of the globe, respectively (Supplementary Table 1). Positive trends in AET and D were jointly observed for 46% of terrestrial surfaces that map onto the PC4 loading. Likewise, joint significant increases in T x,max and T n,min were observed for 24% of the globe and contributed to the positive trend in σ d as T x,max and T n,min were poorly correlated on interannual timescales during the reference period outside of the tropics.
Approximately half of the observed increase in σ d can be accounted for by anthropogenic climate change (Fig. 1e, Supplementary Fig. 2a). The counterfactual simulation produced: 1) σ d trends with a global median of +0.36σ and were spatially similar to that observed, 2) mostly non-significant trends in temperature based departures, and 3) trends in moisture based departures that were approximately a third the magnitude of observed trends. Similarly, excluding the effect of anthropogenic climate change led to a 45% reduction in the extent of global land area with σ d > 2σ (Fig. 1f).
Several factors account for the σ d increase after excluding the modeled influence of anthropogenic climate change, including changes in climate variability and divergence between observed and modeled changes in climate. The sensitivity experiment that excluded trends in annual temperature, precipitation, and reference evapotranspiration had σ d trends that were approximately a third of those observed and trends in climate departures for individual variables that were less than 30% of those observed (Fig. 1e, Supplementary Fig. 1c). By contrast, www.nature.com/scientificreports www.nature.com/scientificreports/ the sensitivity experiment that excluded trends in annual precipitation had nominal influence in attenuating observed trends in climate departures thus highlighting the importance of changes in temperature and evaporative demand (Supplementary Fig. 1b). Our results further suggest that most of the observed increase in climate departures have materialized through changes in mean climate conditions rather than changes in climate variability. Changes in the interannual variability of individual variables showed confounding signals spatially and across variables (Supplemental Methods; Supplementary Table 2). The most notable and widespread change was a global median 15% reduction in the standard deviation of T n,min during 1958-2017. The larger fraction of trends accounted for by explicitly detrending observations versus those using the counterfactual model is expected given that observed σ d trends include changes arising from internal variability.
Lastly, the magnitude of observed changes in climate departures exceeded those estimated from natural variability in a stationary climate. A null model based on 500-years of climate model output run under pre-industrial control simulations 38 yielded a global median σ d trend comparable to that of the detrended sensitivity experiment (average trend +0.24σ over 60-years, maximum of +0.29σ; Supplementary Fig. 6). Internal decadal variability along with the use of a reference period from which climate departures are calculated and associated sampling biases generally results in positive σ d trends even in a stationary climate. However, the substantially larger magnitude of observed changes versus those obtained using the null model suggest that the observed increases in climate departures over our study period are unlikely due to internal climate variability. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
We demonstrate that annual climate conditions over the past three decades have departed substantially from previous decades, that these climate departures exceed what would be expected from natural climate variability, and that approximately half of the magnitude of these changes can be attributed to anthropogenic forcings. While previous studies have highlighted the importance of multivariate extremes and changes thereof 11,13 , we find that multivariate climate departures assessed at annual time scales have rapidly outpaced univariate departures, highlighting the potential latent impacts imposed by observed change to systems adapted to multiple interacting climate variables. The largest positive σ d trends occurred in regions of historically low variance (e.g., equatorial regions) 8,39,40 , in regions such as southern Europe that have seen large changes in climate trends in the variables considered 41 , and in boreal regions that saw joint increases in AET and D that are orthogonal to the historical covariance of these variables. The growing extent of lands with annual multivariate climate departures>2σ in recent decades complements other studies that have considered univariate departures 42 . While our results were specific to the four variables than span moisture and thermal constraints in a given year, trends in multivariate climate departures may differ as a function of the variables chosen.
Trends in multivariate climate departures were largely accounted for by changes in climate means rather than changes in climate variability. While there is a perception that observed climate change has resulted in heightened variability, observational and modeling studies generally do not support widespread detectable changes to date 43,44 . Nonetheless, while observed changes in variability have generally been small, some studies suggest increased hydroclimatic variability or intensification over the 21st century that would further exacerbate increased climate departures [45][46][47] .
Our framing of univariate and multivariate departures builds on previous studies 3,39,48 while emphasizing more nuanced ways in which climate has differed from reference conditions on annual timescales. Only considering changes in mean conditions or changes for multi-decadal timescales can obscure departures that occur annually. This type of interannual variability is critical for assessing climate change impacts. For example, years with concurrent high D and T x,max enable fire activity in flammability-limited forests 19 , while years with high moisture availability (high AET, low D) and low T x,max promote recruitment pulses in water-limited forests 49 . Warm summers accompanied by surface water availability (high T x,max , high AET) may impact human health directly through heat stress and indirectly through promoting environmental conditions that favor certain vector-borne diseases 50,51 . Annual climate departures are also relevant to impacts on agricultural systems, land-use planning, and built infrastructure 52,53 .
As the climate becomes increasingly unfamiliar to the flora and fauna at a specific location, organisms, including humans, must adapt to changing local conditions or move to maintain suitable climates 54 . Indeed, Earth's www.nature.com/scientificreports www.nature.com/scientificreports/ biota is already adapting to climate change. For example, migration and thermophilization of biotic communities has been observed 55,56 . Additionally, humans are increasingly adopting crop varieties, seed provenances, and horticultural practices from other locations with climates that are anticipated to resemble future conditions 3,57 . Unfortunately, the adaptive capacity of systems is not necessarily commensurate with their climate change exposure. For example, thermal specialists in the tropics may be uniquely vulnerable to the large climate departures seen in these regions 58,59 . Collectively, our results suggest that annual multivariate climate departures have changed more dramatically in the last 60 years than univariate estimates of departures for individual climate variables.

Methods
Monthly data from TerraClimate during 1958-2017 were used for our primary analysis 60 . TerraClimate combines high-resolution multi-decadal climatologies from WorldClim 61 with coarser-scale temporally varying data from reanalyses and CRU Ts4.0 62 to create 4-km (1/24°) spatial resolution surfaces of monthly first-order climate variables (e.g., temperature, precipitation), as well as surface water balance products such as AET and D. The water balance model is an accounting-based approach that simulates moisture fluxes using a simple snow model and soil water balance approach for a reference vegetation type. This modeling framework has been widely used in hydrological and ecological studies 10,63 . Other gridded datasets provide similar outputs, including those using more sophisticated surface hydrologic models. However, they do not cover the 60-year period of observational record or are of much coarser spatial resolution. Nonetheless, to interrogate the sensitivity of our results with respect to structural uncertainty of climate data sources, we replicate all analyses for the period 1958-2016 using gridded data from the Princeton Global Meteorological Forcing dataset (version 3) at a 0.25° spatial resolution 64 . Results from these data are presented in Supplemental Methods.
Four climate metrics were chosen given their established links to factors that influence Earth's biota including species occurrence and impacts, and their frequent usage in species distribution models 10,29,32,65 . These metrics included: (1) average maximum temperature of the warmest month for the calendar year (T x,max ), (2) average minimum temperature of the coldest month for the calendar year (T n,min ), (3) calendar year cumulative AET, and (4) calendar year cumulative D. These four variables provide complementary information pertinent to ecological and agricultural systems. For example, T n,min can be a limiting factor for some species and agriculture and modify overwinter morality rates of some organisms 66,67 , while D and AET provide a reduced set of biologically relevant and physically based variables that account for the concurrent availability of both water and energy important for both ecological and agricultural impacts 29,68 .
We calculated σ d and its univariate equivalent standardized Euclidean distances for each variable during 1958-2017. The first 30-years (1958-1987) were used to define a baseline period that determine both the centroid (i.e. the mean) and covariance structure using principal components analysis to calculate σ d (Supplemental Methods). A given 30-year period comprises a sample of a population of climate statistics that may differ from other baseline periods in some regions due to decadal variability 69 . Likewise, anthropogenic climate forcing is evident throughout the observational record limiting the ability to procure a true 'natural' reference state. However, for the basis of this study, the earliest 30-years covers a time period when anthropogenic forcing was substantially less than in recent decades.
We used Mahalanobis distances and its standardization based on the Chi distribution to calculate σ d . This approach accounts for covariance among variables, the dimensionality of the data (number of variables), and local inter-annual climate variability thus facilitating comparisons across space and variable combinations. Mahalanobis distances were calculated on standardized data (i.e., normal distributions based on means and standard deviation calculated during . This approach assumes variables adhere to a multivariate Gaussian distribution, which may be reasonable for some variables, but could be problematic for zero-bounded data (e.g., D) and for temperature extremes in some portions of the globe. We conducted a sensitivity analysis that used clamped non-parametric distributions to assess whether results were substantially altered by these assumptions (Supplemental Methods). The choice of data transformations did not substantially alter the general results of the study, although the truncated nature of the non-parametric distributions reduced overall magnitudes of departures. More advanced multivariate approaches such as copulas may provide additional nuance beyond the linear approaches used herein 70 . Trends for σ d and climate departure for individual variables during 1958-2017 were calculated using Sen-Theil slope estimator and were considered significant using the Mann-Kendall trend test at p < 0.05. Similarly, standard trends using the same approach were calculated on the raw climate variables.
A counterfactual simulation of terrestrial climate was developed that excludes the modeled influence of anthropogenic climate forcing during 1958-2017. The counterfactual simulation removed first-order modeled trends in monthly climate from observed data similar to that of previous studies [71][72][73] . We used a pattern scaling approach that applies the median response among 23 different climate models to proximate modeled anthropogenic changes (Supplemental Methods) 73 . We also considered two sensitivity tests that use detrended observations to decompose changes in climate departures into those associated with trends versus those associated with changes in intraannual-to-interannual variability. One approach removed only annual precipitation trends, while the other approach removed annual trends in temperature, precipitation, and reference evapotranspiration. Detrending was facilitated by calculating a linear least square fit on annual data. This linear fit was subtracted from the observed time series leaving data for 1958 unaltered. By only removing annual trends, we allow monthly trends to differ in a relative sense from annual counterparts. Data for the counterfactual simulation and sensitivity experiments were run through the water balance model to calculate AET and D.
We consider two additional filters in an effort to avoid misinterpreting results in locations of reduced data quality or where the baseline observations were heavily positively-skewed. First, we omit calculations in locations where annual D or AET was 0 for a majority of the years in the baseline period. Second, we used the data quality flags in TerraClimate inherited from CRU TS4.0 to select pixels where at least four stations contributed www.nature.com/scientificreports www.nature.com/scientificreports/ to monthly precipitation, temperature, and vapor pressure fields for at least 75% of the period of record. Fig. S4 shows the data from Fig. 2 with data poor locations masked out. Summarized timeseries and reported trends for global land surfaces exclude pixels not meeting either of these criteria.
Lastly, we use 500-years (model years 400-899) of pre-industrial control climate simulated from the LENS experiment which uses NCAR's Community Earth System Model version 1 (CESM1) with CAM5.2 as its atmospheric model 38 to develop a null model for changes in σ d that result purely from internal climate variability. Monthly output from LENS was subsequently used in the water balance model, albeit at the native spatial resolution of LENS output. We subsequently calculated σ d using non-overlapping moving 30-year blocks over the simulation (e.g., 430-459) and calculated both forward (e.g., 430-489) and backward (e.g., 400-459) linear trends in global terrestrial median σ d and the fraction of land surfaces with σ d > 2σ. The slope of trends calculated from backward samples was inverted. As with observations, trends were calculated using Sen-Theil slope estimator.

Data availability
The primary datasets analyzed in the current study are available through the Northwest Knowledge Network data repository at https://climate.northwestknowledge.net/TERRACLIMATE-DATA/.