Analyzing climate scenarios using dynamic mode decomposition with control

Abstract

Understanding the complex dynamics of climate patterns under different anthropogenic emissions scenarios is crucial for predicting future environmental conditions and formulating sustainable policies. Using Dynamic Mode Decomposition with control (DMDc), we analyze surface air temperature patterns from climate simulations to elucidate the effects of various climate-forcing agents. This improves upon previous DMD-based methods by including forcing information as a control variable. Our study identifies both common climate patterns, like the North Atlantic Oscillation and El Niño Southern Oscillation, and distinct impacts of aerosol and carbon emissions. We show that these emissions’ effects vary with climate scenarios, particularly under conditions of higher radiative forcing. Our findings confirm DMDc’s utility in climate analysis, highlighting its role in extracting modes of variability from surface air temperature while controlling for emissions contributions and exposing trends in these spatial patterns as forcing scenarios change.

Impact Statement

This paper presents an application of dynamic mode decomposition with control to surface air temperature to extract modes of variability and compare different emissions scenarios. It highlights the opportunities of control-based modeling to separately inspect the effect of climate dynamics and forcing from emissions on projections of future climate.

1. Introduction

Modeling and understanding climate modes of variability, like the El Niño-Southern Oscillation (ENSO) and the North Atlantic Oscillation (NAO), are essential for predicting climate impacts on global ecosystems and human activities. Modes of variability can be affected by human activity, significantly influence global weather conditions, and are linked to fundamental physical processes in the atmosphere and ocean. Thus, they are crucial for devising effective climate mitigation and adaptation strategies. For example, modes of variability used to understand the interconnectedness of global weather patterns and their far-reaching effects are often referred to as teleconnections (Diaz et al., 2001; Tsonis et al., 2008; Bridgman and Oliver, 2014).

Dimensionality reduction methods are fundamental in extracting modes of variability for analyzing climate dynamics and teleconnections. Principal Component Analysis (PCA), a.k.a. Empirical Orthogonal Functions (EOF), is traditionally used to extract spatiotemporal features (Bauer-Marschallinger et al., 2013; Volkov, 2014, Forootan et al., 2016). PCA is, however, a linear model, which limits its ability to capture nonlinear processes. Nonlinear formulations, such as kernel PCA, are better suited for nonlinear feature extraction (Schölkopf et al., 1998). Many other PCA extensions have been proposed to address specific climate studies like extended EOF, Multivariate EOF, Principal Oscillation Patterns, and Nonlinear PCA Gehne et al. (2014), Wa et al. (2008), Hannachi et al. (2007), Monahan (2001). Although PCA/EOF and its variants effectively reduce dimensionality, they fall short of capturing the dynamic behavior of climate systems (e.g., the oscillation frequency of modes of variability) due to their inability to process temporal information. This limitation spurred the development of methods like Linear Inverse Modeling (LIM) (Penland and Sardeshmukh, 1995) and Dynamic Mode Decomposition (DMD) (Kutz et al., 2016), which are designed to identify dynamically significant spatiotemporal patterns. LIM continues to find applications in climate science (Wills et al., 2020), though it is not the only dynamical systems-based model used to extract modes of variability in this field. Other methods, such as singular spectrum analysis, are reviewed by Ghil et al. (2002), with subsequent developments including Average Predictability Time (DelSole and Tippett, 2009a, b) and, more recently, Bayesian Linear Dynamical Mode Decomposition (Gavrilov et al., 2020). In contrast, DMD has primarily been applied in fields like fluid dynamics, with a few notable exceptions in climate science, such as emulating sea surface temperature data (Erichson et al., 2019; Navarra et al., 2021), detecting transitions in the North Atlantic Oscillation (Gottwald and Gugole, 2020), and analyzing large-scale climate datasets (Xiong et al., 2023).

Here we use Dynamic Mode Decomposition with control (DMDc) (Proctor et al., 2016), which extends DMD/LIM methods and extracts modes of variability driven by a control variable (Figure 1). Although modes of variability like ENSO indices are usually derived from Sea Surface Temperature (SST) datasets, SSTs undergo trends under climate forcings, which makes it hard to cleanly identify changes in such modes on top of the underlying global warming trend. Our approach applies DMDc to Surface Air Temperature (SAT) data and allows us to separate potential changes in climate modes—such as the ENSO or the NAO—from simultaneous trends in the mean climate state.

Figure 1. The DMDc model for estimating the autoregressive component and emissions contributions to global Surface Air Temperature (SAT). SAT at time $ t $ (in years) is $ \boldsymbol{x} $ , the last $ 30 $ years ( $ t-29,t-28,\dots, t $ ) of radiative forcing are in $ \boldsymbol{y} $ , and SAT at time $ t+1 $ is $ {\boldsymbol{x}}^{\prime } $ . We take $ t=2030 $ in the figure above. The autoregressive component contains the scaled modes of $ \boldsymbol{A} $ , which are scaled eigenvectors of a low-rank estimate of $ \boldsymbol{A} $ . The forcing contribution patterns are entries of a low-rank estimate of $ \boldsymbol{B} $ scaled by entries of $ \boldsymbol{y} $ _.

We highlight the utility of DMDc on SATs under different future emissions scenarios. Specifically, we let $ \boldsymbol{x} $ denote SAT at time $ t $ and $ {\boldsymbol{x}}^{\prime } $ be SAT at time $ t+1 $ . DMD assumes the model $ {\boldsymbol{x}}^{\prime }=\boldsymbol{Ax} $ whereas DMDc includes a linear contribution of the past $ 30 $ years of radiative forcing from emissions in the control variable $ \boldsymbol{y} $ . DMDc results in the model

(1) $$ {\boldsymbol{x}}^{\prime }=\boldsymbol{Ax}+\boldsymbol{By}. $$

This novel application discerns distinct climate responses under multiple emissions scenarios. In particular, DMDc extracts modes of variability that separate into autoregressive modes (eigenvectors of $ \boldsymbol{A} $ ) derived from the year-to-year relationship of SAT with itself and “control” contributions associated with the radiative forcing (using $ \boldsymbol{B} $ ) induced by emissions over the past 30 years. We find that the autoregressive modes include warming trend and ENSO modes, which can be attributed to thermodynamically and dynamically forced trends under different forcing scenarios. We discovered that the emissions contribution from DMDc extracts the variable local impact of radiative forcing for different emissions agents over time, including warming from carbon-based emissions and some cooling from aerosols.

2. Data and methods

This section describes the ClimateBench dataset and summarizes how DMDc is used in our experiments.

2.1. Dataset: ClimateBench

We use the local annual SAT variable from the ClimateBench dataset (Watson-Parris et al., 2022) at approximately 2 degree spatial resolution (96 latitude pixels and 144 longitude pixels), paired with annual emissions for four of the main anthropogenic forcing agents: carbon dioxide (CO₂), methane (CH₄), sulfur dioxide (SO₂) and black carbon (BC). The SAT data is derived from the most recent version of the Norwegian Earth System Model (NorESM2) (Seland et al., 2020) while the emission data is sourced from the Community Emissions Data System (Hoesly et al., 2018).

The data spans the years 1850 to 2050. These years are split into historical simulation from the Coupled Model Intercomparison Project Phase 6 (CMIP6) (Eyring et al., 2016) and different future simulations corresponding to emissions scenarios from ScenarioMIP protocol (O’Neill et al., 2016). The historical data spans the years 1850 to 2014, and the scenario data covers the years 2015 to 2100. We focus on four possible realistic future trajectories called Shared Socioeconomic Pathways (SSPs) (Riahi et al., 2017). From highest to lowest forcing, the selected scenarios are SSP585, SSP370, SSP245, and SSP126 (Figure 2).

Figure 2. The ClimateBench dataset and the derived radiative forcing. This dataset contains four different emissions scenarios: SSP585, SSP370, SSP245, and SSP126 (ordered from highest forcing to lowest forcing). (a) Spatially averaged yearly global temperature from the ClimateBench dataset. (b) Radiative forcing derived from the ClimateBench dataset.

The following analysis uses the SAT from the historical period along with emissions data for each scenario. When these analyses are run on the raw emissions data, forcing agents have a minuscule and inconsistent contribution across emissions scenarios. These contributions do not align with known interactions between forcing agents and surface air temperature. Therefore, we use the procedure outlined in (Leach et al., 2021; Bouabid et al., 2024) to compute estimates of the radiative forcing level induced by each forcing agent from their annual emissions data (Figure 2(b)). Specifically, radiative forcing is computed using the gas cycle model and forcing model of FaIRv2.0.0 calibrated for NorESM2-LM to convert (1) emission data into concentration and (2) concentration into levels of effective radiative forcing (Leach et al., 2021).

2.2. Applying DMDc to ClimateBench

To analyze these data, we use a Python implementation of the DMDc model (Demo et al., 2018; Ichinaga et al., 2024) and specify the radiative forcing contribution from emissions as a control variable. Mirroring the notation by Proctor et al. (2016), we assume the linear model for local SAT in Eq. (1). In this model, $ \boldsymbol{x} $ is the flattened local SAT at year $ t $ , $ \boldsymbol{y} $ is the estimated radiative forcing for each agent from the last 30 years: $ t,t-1,\dots, t-29 $ , and $ {\boldsymbol{x}}^{\prime } $ is the SAT for the next year, $ t+1 $ . Neither the SAT nor the forcing time series are mean-centered before running our analysis. In the DMDc model, $ \boldsymbol{A} $ corresponds to the autoregressive component and $ \boldsymbol{B} $ to the emissions (in this case, radiative forcing) contribution. To simplify the detection of dynamics within the autoregressive component of the system, we assume $ \boldsymbol{A} $ is rank 5. The definitions and symbols for this applying DMDc to ClimateBench appear in Table 1.

Table 1. Definitions of symbols for applying DMDc to SAT and emissions from the ClimateBench dataset

The modes of variability of the autoregressive component, along with their contribution over time, come from the scaled modes, $ {\boldsymbol{\xi}}_m $ , and corresponding eigenvalues, $ {\lambda}_m $ , of $ \boldsymbol{A} $ . The contribution of radiative forcing can be found using the product of entries of $ \boldsymbol{B} $ and $ \mathbf{y} $ corresponding to different forcing agents and times. We now present further details on the DMDc algorithm along with extraction of the autoregressive component and emissions contributions.

2.2.1. Time delayed dataset structure

We first organize the SAT and emissions data from ClimateBench into variables (Table 1).

Our data consists of a flattened local SAT time series

(2) $$ \left\{\boldsymbol{x}(0),\boldsymbol{x}(1),\dots, \boldsymbol{x}(T)\right\}\subset {\mathrm{\mathbb{R}}}^R, $$

where $ R=96\times 144 $ is the dimensionality of the flattened local mean surface temperature grid and a radiative forcing time series (estimated from the ClimateBench emission data following the procedure outlined in (Leach et al., 2021; Bouabid et al., 2024))

(3) $$ \left\{\boldsymbol{y}(0),\boldsymbol{y}(1),\dots, \boldsymbol{y}(T)\right\}\subset {\mathrm{\mathbb{R}}}^S, $$

where $ S=4 $ for the 4 forcing agents considered (CO₂, CH₄, SO₂, and BC).

Let $ \tau =1 $ year be the lag in the autoregressive component (see supplementary material for a detailed explanation), $ J=30 $ years be the delay for the time delay embedding (Takens, 2006) of emissions, and define $ N=T-J-\tau +1 $ . We construct the input and output sample matrices as

(4) $$ \boldsymbol{X}=[\boldsymbol{x}(J)\hskip0.24em \boldsymbol{x}(J+1)\dots \boldsymbol{x}(T-\tau )]\in {\unicode{x211D}}^{R\times N}, $$

(5) $$ {\boldsymbol{X}}^{\prime }=\left[\boldsymbol{x}\left(J+\tau \right)\hskip0.24em \boldsymbol{x}\left(J+\tau +1\right)\dots \boldsymbol{x}(T)\right]\in {\mathrm{\mathbb{R}}}^{R\times N}. $$

This results in a dataset of $ N $ samples, each with $ R $ sample features. Further, using $ K= SJ $ , we define the input control signal following

(6) $$ {\boldsymbol{y}}_n=\left[\begin{array}{c}\boldsymbol{y}(n)\\ {}\boldsymbol{y}\left(n-1\right)\\ {}\vdots \\ {}\boldsymbol{y}\left(n-J+1\right)\end{array}\right]\in {\mathrm{\mathbb{R}}}^K, $$

which we concatenate into the control input matrix

(7) $$ \boldsymbol{Y}=\left[{\boldsymbol{y}}_J\dots {\boldsymbol{y}}_{T-\tau}\right]\in {\mathrm{\mathbb{R}}}^{K\times N}. $$

Thus, we have a dataset of $ N $ control inputs, each with $ K $ control features.

2.2.2. DMDc procedure

Dynamic Mode Decomposition with Control (DMDc) (Proctor et al., 2016) assumes the linear state space model with control

(8) $$ {\boldsymbol{x}}_n^{\prime }={\boldsymbol{Ax}}_n+{\boldsymbol{By}}_n. $$

In matrix form, this model can be block-factorized following

(9) $$ {\boldsymbol{X}}^{\prime }=\boldsymbol{AX}+\boldsymbol{BY}=\left[\boldsymbol{A}\hskip0.24em \boldsymbol{B}\right]\underset{\boldsymbol{\Omega}}{\underbrace{\left[\begin{array}{c}\boldsymbol{X}\\ {}\boldsymbol{Y}\end{array}\right]}}=\left[\boldsymbol{A}\hskip0.24em \boldsymbol{B}\right]\boldsymbol{\Omega} . $$

Using this, DMDc aims to accomplish two tasks.

1. 1. Estimate $ \boldsymbol{A} $ and $ \boldsymbol{B} $ with least-squares regression using the data $ \boldsymbol{X},{\boldsymbol{X}}^{\prime },\boldsymbol{Y} $ .

2. 2. Compute the eigendecomposition of $ \boldsymbol{A} $ , which allows the study of the modes of the autoregressive dynamics.

This is classically achieved using Singular Value Decomposition (SVD) based methods (Tu, 2013), which compute a reduced-rank approximation of $ \boldsymbol{A} $ and $ \boldsymbol{B} $ . The procedure we follow is outlined in Algorithm 1 and outputs eigenvalues along the diagonal of $ \boldsymbol{\Lambda} $ , dynamic modes in the columns of $ \boldsymbol{\Phi} $ , and a reduced rank approximation of $ \boldsymbol{B} $ in . DMDc generalizes DMD by including control (App. A for details).

Algorithm 1 SVD-based DMDc procedure (Proctor et al., 2016).

1: procedure DMDc( $ \boldsymbol{\Omega}, M,{M}_{\Omega} $ )

2: Compute $ {M}_{\Omega} $ -truncated SVD $ \boldsymbol{\Omega} \approx \left[\begin{array}{c}{\boldsymbol{U}}_A\\ {}{\boldsymbol{U}}_B\end{array}\right]\boldsymbol{\Sigma} {\boldsymbol{V}}^{\ast } $

3: Compute $ M $ -truncated SVD $ {\boldsymbol{X}}^{\prime}\approx \hat{\boldsymbol{U}}\hat{\boldsymbol{\Sigma}}{\hat{\boldsymbol{V}}}^{\ast } $

4: Approximate $ \boldsymbol{A} $ as $ \overline{\boldsymbol{A}}={\boldsymbol{X}}^{\prime}\boldsymbol{V}{\boldsymbol{\Sigma}}^{-1}{\boldsymbol{U}}_A^{\ast } $

5: Approximate $ \boldsymbol{B} $ as $ \overline{\boldsymbol{B}}={\boldsymbol{X}}^{\prime}\boldsymbol{V}{\boldsymbol{\Sigma}}^{-1}{\boldsymbol{U}}_B^{\ast } $

6: Compute $ \tilde{\boldsymbol{A}}={\hat{\boldsymbol{U}}}^{\ast}\overline{\boldsymbol{A}}\hat{\boldsymbol{U}} $ and $ \tilde{\boldsymbol{B}}={\hat{\boldsymbol{U}}}^{\ast}\overline{\boldsymbol{B}} $

7: Compute the eigendecomposition $ \left(\boldsymbol{\Lambda}, \boldsymbol{W}\right) $ of $ \tilde{\boldsymbol{A}} $

8: Compute the dynamic modes $ \boldsymbol{\Phi} ={\boldsymbol{X}}^{\prime}\boldsymbol{V}{\boldsymbol{\Sigma}}^{-1}{\boldsymbol{U}}_A^{\ast}\hat{\boldsymbol{U}}\boldsymbol{W} $

9: Rank- $ M $ approximation

10: return .

11: end procedure

In experiments, we use a full-rank SVD of $ \boldsymbol{\Omega} $ , e.g., $ {M}_{\Omega}=\operatorname{rank}\left(\boldsymbol{\Omega} \right) $ . We compute a reduced-rank approximation of $ \boldsymbol{B} $ by restricting the rank of the SVD of $ {\boldsymbol{X}}^{\prime } $ . This chosen rank $ M $ will determine the number of dynamic modes. We choose $ M=5 $ as we experimentally observe it captures the dominant singular values (see supplementary material). Instead of using the least squares approximations ( $ \overline{\boldsymbol{A}} $ and $ \overline{\boldsymbol{B}} $ ) to emulate or reconstruct the SAT series, we use the first $ M=5 $ dynamic modes from Algorithm 1 and to extract modes of variability.

2.2.3. Autoregressive component

We use DMDc to analyze the dynamics from the matrix $ \boldsymbol{A} $ by first computing the eigenvalue decomposition $ \tilde{\boldsymbol{A}}\boldsymbol{W}=\boldsymbol{W}\boldsymbol{\Lambda } $ . This eigenvalue decomposition gives rise to the eigenvalues, dynamic modes, and amplitudes (Table 2). The spatial pattern is captured in the dynamic mode, the amplitude captures the fixed contribution of the dynamic mode, and the eigenvalue summarizes the oscillation and trend over time.

Table 2. Decomposition of $ \tilde{\boldsymbol{A}} $

Spatial patterns. The vector of amplitudes, $ \boldsymbol{\alpha} ={[{\alpha}_1\cdots {\alpha}_M]}^T $ , is defined as the minimizer of the least squares problem $ \boldsymbol{\Phi} \boldsymbol{\alpha} \approx {\boldsymbol{x}}_1 $ . The product of a dynamic mode and its corresponding, fixed amplitude is the scaled mode (Krake et al., 2021) denoted

(10) $$ {\boldsymbol{\xi}}_m={\alpha}_m{\boldsymbol{\varphi}}_m. $$

We visualize the scaled modes in our experiments to provide a realistic scaling (e.g., correcting for sign flips) of the spatial patterns of the dynamic modes (Figures 4, 5, and 6).

Temporal patterns. Inspired by (Krake et al., 2021), we determine the evolution of the scaled modes over time with the help of the following propositions. See supplementary material for the proofs.

Proposition 1. The dynamic modes, $ {\boldsymbol{\varphi}}_m $ , are the eigenvectors of $ {\overline{\boldsymbol{A}}\hat{\boldsymbol{U}}\hat{\boldsymbol{U}}}^{\ast } $ with eigenvalues $ {\lambda}_m $ .

Proposition 2. $ {\lambda}_m{\boldsymbol{\xi}}_m $ is the evolution of a scaled mode, $ {\boldsymbol{\xi}}_m $ , over $ \tau $ time steps in the subspace spanned by the columns of $ \hat{\boldsymbol{U}} $ . Specifically,

(11) $$ {\lambda}_m{\boldsymbol{\xi}}_m={\overline{\boldsymbol{A}}\hat{\boldsymbol{U}}\hat{\boldsymbol{U}}}^{\ast }{\boldsymbol{\xi}}_m\in {\mathrm{\mathbb{R}}}^R. $$

By linearity, the evolution of two scaled modes over $ \tau $ time steps is $ {\lambda}_j{\boldsymbol{\xi}}_j+{\lambda}_k{\boldsymbol{\xi}}_k. $ We use this fact to visualize the evolution of two scaled modes with complex conjugate eigenvalues (Figure 6b).

The trend and oscillation frequency of a scaled mode can be extracted from its associated eigenvalue. When written in exponential form, the $ m $ th eigenvalue is $ {\lambda}_m={r}_m{e}^{i{\omega}_m} $ . The magnitude $ {r}_m=\mid {\lambda}_m\mid $ is the magnitude of the contribution of the scaled mode over time. A magnitude of $ {r}_m=1 $ is stable, $ {r}_m<1 $ is decaying, and $ {r}_m>1 $ is a diverging contribution. The complex portion $ {e}^{i{\omega}_m} $ controls the frequency of oscillation

(12) $$ {f}_m=\frac{\omega_m}{2\pi }=-i\log \left(\frac{\lambda_m}{\mid {\lambda}_m\mid}\right)/\left(2\pi \right). $$

We combine the previous two propositions and the definition of the dynamic mode amplitudes. Using $ \boldsymbol{\Xi} $ to represent the matrix whose columns are scaled dynamic modes ( $ {\boldsymbol{\xi}}_m $ ), this results in the approximation of the autoregressive part as.

$ {\boldsymbol{Ax}}_n\approx \boldsymbol{\Xi} \Lambda {\boldsymbol{\Xi}}^{\dagger }{\boldsymbol{x}}_n. $

The details of this approximation can be found in the supplementary material.

2.2.4. Emissions contribution

The columns of the matrix are fixed spatial patterns of the linear forcing contribution from emissions (Figure 3). The scale of each of these contributions to SAT is exactly the corresponding entry in $ \boldsymbol{y} $ .

Figure 3. The structure of the forcing contribution matrix $ \overset{\smile }{\boldsymbol{B}} $ _.

The evolution of the contribution of radiative forcing over times $ \left(1,2,\dots, N\right) $ is . The mean of each vector is the spatial mean contribution of emissions to temperature over time (Figure 7 top left).

Time-lagged radiative forcing. We use the time-lagged structure of $ {\boldsymbol{y}}_n $ to separate the radiative forcing signal into years (Eq. (6)). Each entry of $ \boldsymbol{y}(j)\in {\mathrm{\mathbb{R}}}^S $ (denoted $ \boldsymbol{y}{(j)}_s $ ) corresponds to a different forcing agent (e.g., CO₂). In our experiments, we have $ S=4 $ forcing agents.

Structure of _. The spatial contribution of radiative forcing over the last $ J $ years is determined by decomposing using the structure of $ {\boldsymbol{y}}_n $ . First, is decomposed into blocks of size $ N\times S $ corresponding to each year in $ {\boldsymbol{y}}_n $

(13)

Then the spatial contribution of the input radiative forcing agent $ s $ (e.g., CO₂) over $ \left(1,2,\dots, J\right) $ time steps in the past are grouped into the $ {s}^{\mathrm{th}} $ columns of each of , and is denoted . This structure is summarized in Figure 3.

This process unravels the matrix multiplication into pieces corresponding to time lag and forcing agent as

(14)

The spatial pattern corresponding to the linear contribution of a forcing agent $ s $ at year $ j $ before the predicted year $ n+1 $ is found in the entries of . Using this method, we determine the forcing contribution for each agent as the products between entries of and $ {\boldsymbol{y}}_n $ (Figures 7–9). See Appendix B for details of the usage.

Limitations. DMDc does not explicitly address the possible correlation between the time series dimensions and, as a result, can suffer from the same pitfalls as linear regression (e.g., difficulty in interpreting the coefficients because they may not linearly separate). If we are using a form of regularization in fitting the $ \boldsymbol{A} $ and $ \boldsymbol{B} $ matrices, that may help mitigate this effect by ensuring we do not “mix up” factors as much. We leave this as an interesting avenue for future study. We acknowledge that this potential for correlation limits the extent to which we state that one particular factor is responsible for the obtained pattern. For example, a correlation between CO₂ and CH₄ forcings can be due to a shared latent variable (human activity) that drives them, and interpreting the coefficients of a linear model such as DMDc as the causal effect of a forcing (CO₂ or CH₄) has on the pattern does not correct for this confounding.

Summary. Our application of DMDc to SAT results in the decomposition of SAT at year $ n+1 $ as

(15)

In this decomposition, the $ {m}^{\mathrm{th}} $ column of $ \boldsymbol{\Xi} $ is a scaled mode ( $ {\boldsymbol{\xi}}_m $ ) and the $ {m}^{\mathrm{th}} $ diagonal entry of $ \boldsymbol{\Lambda} $ is its associated eigenvalue ( $ {\lambda}_m $ ). The term is the spatial pattern of the forcing contribution of emissions agent $ s $ at $ j $ years before year $ n+1 $ .

3. Results

For our analysis, we use DMDc on the ClimateBench dataset (Watson-Parris et al., 2022) to detect modes of variability from local, gridded annual mean SAT at $ \approx 2{}^{\circ} $ horizontal resolution while including the forcing information from various emissions modalities (Figure 1 and Section 2). We compare these representations across four emissions scenarios: the Shared Socio-economic Pathways (SSPs) 585, 370, 245, and 126. This analysis results in modes of variability from (1) the autoregressive component: year-to-year autoregression of the SAT time series (from $ \boldsymbol{A} $ ), and (2) the emissions contribution: the effect of the last 30 years of radiative forcing from different emissions types on SAT (from $ \boldsymbol{B} $ ).

3.1. Autoregressive component

The scaled modes of the autoregressive component of DMDc detect the year-to-year autoregressive contribution of SAT. We use scaled modes to remove sign ambiguity for spatial plots (Section 2). This component contains the effect of the previous year’s temperature on the next year and the forcing contributions that are not captured by emissions contribution. Scaled modes are fixed spatial patterns that oscillate over time and are extracted from the data fields using the eigenvectors of a low-rank estimate of $ \boldsymbol{A} $ . The temporal evolution of a scaled mode is determined by its associated eigenvalue (Section 2).

The scaled modes extracted from SAT display a clear warming trend and other common climate patterns, effectively capturing the year-to-year impact of previous temperatures and external forcings (Figures 4–6) (Hasselmann, 1976). DMD run on the same data extracts similar spatial patterns, trends, and oscillations (see comparison to DMD in supplementary material). However, the DMD model does not incorporate emissions information and thus cannot extract the essential emissions contribution patterns detected by DMDc.

Figure 4. The global warming trend is associated with the two largest real eigenvalues. The scenarios ordered from high to low emissions are SSP585 (blue), SSP370 (orange), SSP245 (green), and SSP126 (red). Each spatial pattern is the sum of the two scaled dynamic modes associated with the two largest real eigenvalues. The plot on the right is the mean latitudinal temperature profile. The color bars and latitudinal profiles are in °C.

The global warming trend is captured in two scaled modes related to warming and cooling with associated eigenvalues. We plot the sum of these two scaled modes in Figure 4. Higher emissions scenarios have a higher magnitude warming trend. Higher emission scenarios have larger eigenvalues for the warming mode (0.98 for SSP585 and 0.93 for SSP126), indicating increased stability over time.

Global warming modes reveal spatial patterns such as intensified land warming in the Northern Hemisphere, especially in Russia (Figure 5), and the North Atlantic Warming Hole (NAWH) (Collins et al., 2013) (Figure 4). There are various mechanisms under discussion to contribute to the existence of the NAWH, from a weakened Atlantic Meridional Overturning Circulation (AMOC) to changes in atmospheric wind patterns (Drijfhout et al., 2012; Bellomo et al., 2021; He et al., 2022) as a combination of the effects of aerosol and greenhouse gas forcing (Qasmi, 2023). North Atlantic temperature variability is generally also associated with the North Atlantic Oscillation (NAO) due to two-way interactions between the atmosphere and ocean (Wang et al., 2004). Since the warming pattern does not oscillate, this suggests that the NAWH can be interpreted as a thermodynamically forced trend rather than a feature of long-term internal variability in the NAO. Additionally, all scenarios feature a South Atlantic warming hole that exhibits no discernible trend with rising emissions. Moreover, all scenarios show a mild El Niño-like signal. This pattern is often connected to a weakening of the Pacific zonal Walker circulation under CO₂ forcing in climate models, which enforces our discovery that this pattern is weakest in the high aerosol forcing scenario, SSP370 (Watanabe et al., 2024).

Figure 5. The SAT warming trend by land mass region over the entire scenario. The scenarios ordered from high to low emissions are SSP585, SSP370, SSP245, and SSP126. Each boxplot summarizes the distribution of the global warming scaled mode values in each region. The whiskers cover the entire distribution of the temperature values. We see larger warming relative to the global mean for the higher forcing scenarios in the Northern Hemisphere, especially Russia.

Two other scaled modes are associated with complex eigenvalues indicating oscillation over time (Figure 6). The tropical Pacific patterns we find in these modes have their peak positive amplitude in the Central Pacific—somewhat akin to an average Central Pacific El Niño event. While our approach to characterizing climate mode behavior is different, our pattern resembles those found in the ENSO analysis in the original NorESM2 model documentation paper by Seland et al., (2020, Figure 23). There, the El Niño composite signal of NorESM2 showed a larger amplitude than observed, with a clear peak amplitude in the Central Pacific (rather than a more elongated East-to-Central Pacific warm anomaly). Therefore, we call these scaled modes the “ENSO modes.” These modes oscillate and have a clear pattern that resembles central Pacific El Niño events for all scenarios except SSP585, which is a La Niña pattern that evolves into a central Pacific El Niño event one year into the future using the associated eigenvalues. These ENSO patterns oscillate at somewhat different frequencies for different scenarios (from 0.15/year for SSP126 to 0.18/year for SSP585). Except for SSP245, the frequency of ENSO increases in higher emissions scenarios (see Autoregressive component in supplementary material). There is evidence of this pattern of increased emissions in conjunction with increasing ENSO variability in many climate model ensembles (Timmermann et al., 1999; Malik et al., 2020). However, these oscillation differences could be estimation uncertainty given the application of DMDc to single ensemble members.

Figure 6. Analysis of ENSO oscillation with DMDc. The scenarios ordered from high to low emissions are SSP585(blue), SSP370 (orange), SSP245 (green), and SSP126 (red). The plot on the right of each map is the mean latitudinal temperature profile. The color bars and latitudinal profiles are in °C. (a) ENSO-like spatial patterns appear in all emissions scenarios, specifically, a central Pacific El Niño event in SSP126, SSP245, and SSP370, whereas La Niña in SSP585. Each spatial pattern is the sum of the two scaled dynamic modes associated with complex conjugate eigenvalues. (b) We observe a Niña-to-Niño phase transition in SSP585 as we evolve the spatial pattern one year into the future using the associated eigenvalues (Section 2). Therefore, a pattern akin to a central Pacific El Niño event is visible in all scenarios. (Left) The spatial pattern is the sum of the two scaled dynamic modes associated with complex conjugate eigenvalues. (Right) The spatial pattern is the sum of the two scaled dynamic modes associated with complex conjugate eigenvalues, multiplied by their associated eigenvalue.

3.2. Emissions contribution

The emissions contribution contains the linear impact of the past 30 years of radiative forcing from the four of the main anthropogenic forcing agents: carbon dioxide (CO₂), methane (CH₄), sulfur dioxide (SO₂) and black carbon (BC) on the current year SAT. Recall that these patterns cannot be detected using methods like DMD because they do not include emissions information. The emissions contribution includes part of the linear thermodynamic response of SAT to different emissions modalities because, in DMDc, they linearly shift the next year’s temperature.

The spatial mean of the emissions contribution increases for higher emissions scenarios (Figure 7, top left). The impact of each of the last 30 years of forcing on a current year increases as we move further into the past (Figure 7, top right). Carbon radiative forcing dominates the spread of the forcing contribution (Figure 7, bottom row) and is known to have a delayed effect on temperature (Zickfeld and Herrington, 2015). This explains the increased impact of radiative forcing on SAT as we look further into the past. Additionally, we compute the grid point standard deviation over all forcing agents from 2051 to 2100. We see an increase in standard deviation from 0.14 °C in SSP126 to 0.20 °C in SSP585 which indicates an overall higher variability in forcing impact in higher forcing scenarios.

Figure 7. Analysis of emissions contribution with DMDc. The scenarios ordered from high to low emissions are SSP585, SSP370, SSP245, and SSP126. Units are °C for the vertical axis in all plots. (Top left) Each annual value represents the spatial mean of the past 30 years of radiative forcing contributions from emissions to global SAT. (Top right) The spatial mean of the cumulative forcing from years 2051–2100 as we look further back from the predicted year averaged over across all four scenarios. (Bottom) The distribution of SAT contribution for each forcing agent for each scenario. The whiskers include each point in the data distribution.

Although DMDc properly isolates the effects of radiative forcing from carbon on SAT, it does not fully capture a mean global cooling from aerosols (see forcing contribution in supplementary material). However, the forcing contribution does capture cooling from SO₂ on northern landmasses with high human populations under some scenarios (SSP585 in Figure 8). The warming in other scenarios could be driven by relative change due to decreasing SO₂ emissions over time (Figure 2(b) in Section 2 for the decreasing magnitude of SO₂).

Figure 8. The forcing contribution to SAT from SO₂ by land mass region looking 30 years backwards from the year 2101. Each boxplot is the distribution of SAT contribution for the forcing contribution of SO₂. The whiskers cover each point in the data distribution. The scenarios, ordered from high to low emissions, are SSP585, SSP370 SSP245 and SSP126 Although SSP585 exhibits cooling from SO₂, we see warming from SO₂ in other scenarios.

The local impacts of the total forcing contribution include increased warming in high latitudes of the Northern Hemisphere (Figure 9). This polar warming pattern is consistent with Arctic amplification, a well-documented phenomenon in climate science literature: polar regions warm faster than the rest due principally to ice-albedo feedback mechanisms (Pithan and Mauritsen, 2014; Rantanen et al., 2022). This pattern is especially expected under the high greenhouse gas emissions (Meehl et al., 2020).

Figure 9. Visualization of the linear impact of radiative forcing on SAT. The scenarios ordered from high to low emissions are SSP585 (blue), SSP370 (orange), SSP245 (green), and SSP126 (red). The plot on the right of each map is the mean latitudinal temperature profile. The color bars and latitudinal profiles are in °C. (a) The cumulative effect of the past 30 years of radiative forcing on SAT (averaged over the output years 2050 to 2100). We see an increased contribution to SAT for higher forcing scenarios and polar amplification. Due to the DMDc model, this only contains the linear contribution of radiative forcing on SAT. (b) The changing pattern of the effect of radiative forcing from carbon on SAT for different years before the predicted year in SSP585 averaged over predicted years $ 2050 $ to $ 2100 $ ). We see a higher contribution from emissions when we look further into the past from the predicted year.

4. Conclusion

We used DMDc to extract modes of variability and emissions contributions under different climate scenarios. To do this, we leveraged ClimateBench, a novel, rich dataset that provides local, gridded yearly mean surface air temperatures while including crucial emissions information. Previous methods (e.g., LIM/DMD) are limited to only considering an autoregressive model for surface air temperature while ignoring essential control information like radiative forcing. DMDc addresses this issue by incorporating radiative forcing as a control variable. When applied to surface air temperatures, DMDc determined distinct (e.g., dynamic versus thermodynamic) climate responses under different emissions scenarios, producing insights into climate change dynamics.

Spatially explicit patterns and temporally resolved trends of impacts were quantified through running DMDc on these data. Specifically, we detected a global warming trend and ENSO (an important modulator of climate extremes) in the autoregressive component and also quantified their increase in higher emissions scenarios. We uncovered the known North Atlantic warming hole and saw its strength decrease as we increased the emissions scenario. Unique to DMDc, the forcing component extracted spatial contributions of four radiative forcing agents: CO₂, CH₄, SO₂, and BC over the last 30 years. As expected, the forcing contribution is higher in higher emissions scenarios. We also saw warming from carbon across all scenarios and some landmass cooling from aerosol in specific scenarios. Other, more nuanced patterns identified by DMDc include changes in emissions contributions as we analyze the effect of emissions over the past 30 years. Finally, we identified spatial patterns of the forcing contributions, which included Arctic amplification.

Our future work will include algorithm refinement and application of DMDc to other climatological variables. For example, our method for determining the emissions contribution from DMDc does not account for potential partial correlations between variables. Exploration of non-linear forcing contribution and other DMD variants (e.g., adding constraints) may address partial correlations, increase forcing contribution to SAT (Deng et al., 2020), improve extraction of the aerosol-based cooling signal, and result in more temporally stable (e.g., non-decaying) autoregressive modes. Finally, further analysis of the control signal with larger datasets could lead to a novel detection attribution strategy. Overall, our study demonstrates the opportunities of DMDc for climate science, paving the way for more novel, data-driven analysis of climate data.