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., Reference Diaz, Hoerling and Eischeid2001; Tsonis et al., Reference Tsonis, Swanson and Wang2008; Bridgman and Oliver, Reference Bridgman and Oliver2014).
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., Reference Bauer-Marschallinger, Dorigo, Wagner and Van Dijk2013; Volkov, Reference Volkov2014, Forootan et al., Reference Forootan, Khandu, Awange, Schumacher, Anyah, Van Dijk and Kusche2016). 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., Reference Schölkopf, Smola and Müller1998). 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. (Reference Gehne, Kleeman and Trenberth2014), Wa et al. (Reference Wa, Wu, Li, Liu, Chang, Ding and Wu2008), Hannachi et al. (Reference Hannachi, Jolliffe and Stephenson2007), Monahan (Reference Monahan2001). 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, Reference Penland and Sardeshmukh1995) and Dynamic Mode Decomposition (DMD) (Kutz et al., Reference Kutz, Brunton, Brunton and Proctor2016), which are designed to identify dynamically significant spatiotemporal patterns. LIM continues to find applications in climate science (Wills et al., Reference Wills, Sippel and Barnes2020), 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. (Reference Ghil, Allen, Dettinger, Ide, Kondrashov, Mann, Robertson, Saunders, Tian and Varadi2002), with subsequent developments including Average Predictability Time (DelSole and Tippett, Reference DelSole and Tippett2009a, Reference DelSole and Tippettb) and, more recently, Bayesian Linear Dynamical Mode Decomposition (Gavrilov et al., Reference Gavrilov, Kravtsov and Mukhin2020). 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., Reference Erichson, Mathelin, Kutz and Brunton2019; Navarra et al., Reference Navarra, Tribbia and Klus2021), detecting transitions in the North Atlantic Oscillation (Gottwald and Gugole, Reference Gottwald and Gugole2020), and analyzing large-scale climate datasets (Xiong et al., Reference Xiong, Ma, Sun and Tian2023).
Here we use Dynamic Mode Decomposition with control (DMDc) (Proctor et al., Reference Proctor, Brunton and Kutz2016), 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
$$ {\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., Reference Watson-Parris, Rao, Olivié, Seland, Nowack, Camps-Valls, Stier, Bouabid, Dewey and Fons2022) 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 (CO2), methane (CH4), sulfur dioxide (SO2) and black carbon (BC). The SAT data is derived from the most recent version of the Norwegian Earth System Model (NorESM2) (Seland et al., Reference Seland, Bentsen, Olivié, Toniazzo, Gjermundsen, Graff, Debernard, Gupta, He and Kirkevåg2020) while the emission data is sourced from the Community Emissions Data System (Hoesly et al., Reference Hoesly, Smith, Feng, Klimont, Janssens-Maenhout, Pitkanen, Seibert, Vu, Andres and Bolt2018).
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., Reference Eyring, Bony, Meehl, Senior, Stevens, Stouffer and Taylor2016) and different future simulations corresponding to emissions scenarios from ScenarioMIP protocol (O’Neill et al., Reference O’Neill, Tebaldi, Van Vuuren, Eyring, Friedlingstein, Hurtt, Knutti, Kriegler, Lamarque and Lowe2016). 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., Reference Riahi, Van Vuuren, Kriegler, Edmonds, O’neill, Fujimori, Bauer, Calvin, Dellink and Fricko2017). 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., Reference Leach, Jenkins, Nicholls, Smith, Lynch, Cain, Walsh, Wu, Tsutsui and Allen2021; Bouabid et al., Reference Bouabid, Sejdinovic and Watson-Parris2024) 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., Reference Leach, Jenkins, Nicholls, Smith, Lynch, Cain, Walsh, Wu, Tsutsui and Allen2021).
2.2. Applying DMDc to ClimateBench
To analyze these data, we use a Python implementation of the DMDc model (Demo et al., Reference Demo, Tezzele and Rozza2018; Ichinaga et al., Reference Ichinaga, Andreuzzi, Demo, Tezzele, Lapo, Rozza, Brunton and Kutz2024) and specify the radiative forcing contribution from emissions as a control variable. Mirroring the notation by Proctor et al. (Reference Proctor, Brunton and Kutz2016), 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
$$ \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., Reference Leach, Jenkins, Nicholls, Smith, Lynch, Cain, Walsh, Wu, Tsutsui and Allen2021; Bouabid et al., Reference Bouabid, Sejdinovic and Watson-Parris2024))
$$ \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 (CO2, CH4, SO2, 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, Reference Takens2006) of emissions, and define
$ N=T-J-\tau +1 $
. We construct the input and output sample matrices as
$$ \boldsymbol{X}=[\boldsymbol{x}(J)\hskip0.24em \boldsymbol{x}(J+1)\dots \boldsymbol{x}(T-\tau )]\in {\unicode{x211D}}^{R\times N}, $$
$$ {\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
$$ {\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
$$ \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., Reference Proctor, Brunton and Kutz2016) assumes the linear state space model with control
$$ {\boldsymbol{x}}_n^{\prime }={\boldsymbol{Ax}}_n+{\boldsymbol{By}}_n. $$
In matrix form, this model can be block-factorized following
$$ {\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. Estimate
$ \boldsymbol{A} $
and
$ \boldsymbol{B} $
with least-squares regression using the data
$ \boldsymbol{X},{\boldsymbol{X}}^{\prime },\boldsymbol{Y} $
. -
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, Reference Tu2013), 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., Reference Proctor, Brunton and Kutz2016).
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., Reference Krake, Reinhardt, Hlawatsch, Eberhardt and Weiskopf2021) denoted
$$ {\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., Reference Krake, Reinhardt, Hlawatsch, Eberhardt and Weiskopf2021), 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,
$$ {\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
$$ {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., CO2). 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 $
Then the spatial contribution of the input radiative forcing agent
$ s $
(e.g., CO2) 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
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 CO2 and CH4 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 (CO2 or CH4) 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
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., Reference Watson-Parris, Rao, Olivié, Seland, Nowack, Camps-Valls, Stier, Bouabid, Dewey and Fons2022) 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, Reference Hasselmann1976). 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., Reference Collins, Knutti, Arblaster, Dufresne, Fichefet, Friedlingstein, Gao, Gutowski, Johns and Krinner2013) (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., Reference Drijfhout, Van Oldenborgh and Cimatoribus2012; Bellomo et al., Reference Bellomo, Angeloni, Corti and von Hardenberg2021; He et al., Reference He, Clement, Cane, Murphy, Klavans and Fenske2022) as a combination of the effects of aerosol and greenhouse gas forcing (Qasmi, Reference Qasmi2023). 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., Reference Wang, Anderson, Kaufmann and Myneni2004). 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 CO2 forcing in climate models, which enforces our discovery that this pattern is weakest in the high aerosol forcing scenario, SSP370 (Watanabe et al., Reference Watanabe, Kang, Collins, Hwang, McGregor and Stuecker2024).
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., (Reference Seland, Bentsen, Olivié, Toniazzo, Gjermundsen, Graff, Debernard, Gupta, He and Kirkevåg2020, 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., Reference Timmermann, Oberhuber, Bacher, Esch, Latif and Roeckner1999; Malik et al., Reference Malik, Nowack, Haigh, Cao, Atique and Plancherel2020). 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 (CO2), methane (CH4), sulfur dioxide (SO2) 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, Reference Zickfeld and Herrington2015). 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 SO2 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 SO2 emissions over time (Figure 2(b) in Section 2 for the decreasing magnitude of SO2).
Figure 8. The forcing contribution to SAT from SO2 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 SO2. 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 SO2, we see warming from SO2 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, Reference Pithan and Mauritsen2014; Rantanen et al., Reference Rantanen, Karpechko, Lipponen, Nordling, Hyvärinen, Ruosteenoja, Vihma and Laaksonen2022). This pattern is especially expected under the high greenhouse gas emissions (Meehl et al., Reference Meehl, Senior, Eyring, Flato, Lamarque, Stouffer, Taylor and Schlund2020).
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: CO2, CH4, SO2, 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., Reference Deng, Dai and Xu2020), 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.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/eds.2025.8.
Acknowledgements
N.M. thanks H. Durand for illuminating discussions about DMDc and climate science.
Author contribution
N.M. and S.B. designed the study with the help of G.C-V. N.M. performed the study and wrote the manuscript with the help of all coauthors. G.C-V. supervised the study. All authors reviewed the manuscript.
Competing interest
The authors declare no competing interests.
Data availability statement
The data used in this manuscript can be found here https://zenodo.org/records/7064308 under version 1.0.0 in the files test.tar.gz and train_val.tar.gz. The specific files used are of the form inputs_ssp245.nc, outputs_ssp245.nc, inputs_historical.nc, and outputs_historical.nc. Code to reproduce the results can be found here https://zenodo.org/records/14730900.
Ethical standard
The research meets all ethical guidelines, including adherence to the legal requirements of the study country.
Funding statement
N.M., D.B., and G.C-V. acknowledge the support of Generalitat Valenciana and the Conselleria d’Innovació, Universitats, Ciéncia i Societat Digital, through the project “AI4CS: Artificial Intelligence for complex systems: Brain, Earth, Climate, Society” (CIPROM/2021/56). G.C-V. acknowledges the support from the European Research Council (ERC) under the ERC Synergy Grant USMILE (grant agreement 855,187) and the HORIZON program under the AI4PEX project (grant agreement 101,137,682). P.N. was supported by the UK Natural Environment Research Council (grant no. NE/V012045/1).
A. From DMD/LIM to DMDc
In DMD (Tu, Reference Tu2013) (a.k.a. LIM (Penland and Sardeshmukh, Reference Penland and Sardeshmukh1995)), the model is
$ d\boldsymbol{x}/ dt=\boldsymbol{Cx}+\boldsymbol{\varepsilon} $
where
$ \boldsymbol{\varepsilon} $
is noise. Without the noise term and a time lag
$ \tau =1 $
, this translates to
$$ \boldsymbol{x}(t+1)=\boldsymbol{A}\boldsymbol{x}(t)=\mathrm{exp}(\boldsymbol{C})\boldsymbol{x}(t). $$
Generally, LIMs use this model where the features of
$ \boldsymbol{x}(t) $
are the principal components of the data rather than the raw data. LIMs are used as naive estimators of the forced response from climate variables via the “least damped mode” of
$ \boldsymbol{A} $
(Solomon and Newman, Reference Solomon and Newman2012). This corresponds to the mode of
$ \boldsymbol{A} $
associated with the largest, entirely real eigenvalue. If the forced response is in the least damped mode, the internal variability is thus contained in the other modes of
$ \boldsymbol{A} $
.
We will now re-capitulate (projected) DMD (Tu, Reference Tu2013). Define the snapshot matrices as before, and take the rank
$ M $
SVD:
$ \boldsymbol{X}\approx \boldsymbol{U}\boldsymbol{\Sigma } {\boldsymbol{V}}^{\ast } $
. DMD assumes that the data satisfy
$ {\boldsymbol{U}}^{\ast }{\boldsymbol{x}}_{n+1}=\tilde{\boldsymbol{A}}{\boldsymbol{U}}^{\ast }{\boldsymbol{x}}_n $
where
$ \tilde{\boldsymbol{A}}={\boldsymbol{U}}^{\ast}\tilde{\boldsymbol{A}}={\boldsymbol{U}}^{\ast }{\boldsymbol{X}}^{\prime}\boldsymbol{V}{\boldsymbol{\Sigma}}^{-1} $
with eigenvectors
$ \tilde{\boldsymbol{A}}\boldsymbol{W}=\boldsymbol{W}\boldsymbol{\Lambda } $
. We use the (projected) dynamic modes,
$ \boldsymbol{\Phi} =\boldsymbol{UW} $
. Define the amplitudes as the least squares approximation
$ \boldsymbol{\Phi} \boldsymbol{\alpha} \approx {\boldsymbol{x}}_1 $
and the scaled modes as
$ {\boldsymbol{\xi}}_m={\boldsymbol{\varphi}}_m{\alpha}_m $
. As with DMDc, the scaled modes make up the autoregressive contribution. Multiplication of a scaled mode by its associated eigenvalue,
$ {\lambda}_m $
evolves the scaled mode one-time step into the future. The DMDc model generalizes LIM/DMD by including a linear control term instead of noise. The DMDc model with a time lag of
$ \tau =1 $
is
$ {\hat{\boldsymbol{U}}}^{\ast }{\boldsymbol{x}}_{n+1}={\tilde{\boldsymbol{A}}\hat{\boldsymbol{U}}}^{\ast }{\boldsymbol{x}}_n+\tilde{\boldsymbol{B}}{\boldsymbol{y}}_n $
. A summary of these differences is in Table A1.
Table A1. The spatial and temporal patterns from DMD and DMDc for climate analysis. Both methods have spatial and temporal patterns for the autoregressive component (scaled modes of
$ \boldsymbol{A} $
). DMDc adds spatial and temporal patterns for the emissions contribution by analyzing the matrix
(Figure 3) using the time-lagged radiative forcing stored in
$ \boldsymbol{y} $
B. Emissions contribution
In the upper left panel of Figure 7, we plot the mean entry of
. Let
$ \mathcal{N} $
be the set of all
$ n $
corresponding to years 2051 to 2100. In the upper right panel, we plot the spatial mean entry of
for
$ j=1,2,\dots, 30 $
years in the past and emissions agent
$ s $
. In the bottom panel, we visualize the distribution of the entries of
over
$ j=1,2,\dots, 30 $
years in the past and a fixed emissions agent
$ s $
. In Figure 8, we plot the distribution of the entries of
for
$ s $
corresponding to SO2 and certain landmass regions. In Figure 9(a) we plot
In Figure 9(b) we plot
for
$ j=1 $
and
$ j=30 $
.
Open access
