# Detecting changepoints
in globally-indexed functional time series

Drew Yarger and J. Derek Tucker
Statistical Sciences, Sandia National Laboratories

###### Abstract

In environmental and climate data, there is often an interest in determining
if and when changes occur in a system. Such changes may result from localized
sources in space and time like a volcanic eruption or climate geoengineering
events. Detecting such events and their subsequent influence on climate has
important policy implications. However, the climate system is complex, and
such changes can be challenging to detect. One statistical perspective for
changepoint detection is functional time series, where one observes an entire
function at each time point. We will consider the context where each time
point is a year, and we observe a function of temperature indexed by day of
the year. Furthermore, such data is measured at many spatial locations on
Earth, which motivates accommodating sets of functional time series that are
spatially-indexed on a sphere. Simultaneously inferring changes that can occur
at different times for different locations is challenging. We propose test
statistics for detecting these changepoints, and we evaluate performance using
varying levels of data complexity, including a simulation study, simplified
climate model simulations, and climate reanalysis data. We evaluate changes in
stratospheric temperature globally over 1984-1998. Such changes may be
associated with the eruption of Mt. Pinatubo in 1991.

Keywords: atmospheric science, changepoint analysis, spatially-dependent
functional data

##  1 Introduction

In environmental and climate data, there is often an interest in determining
if and when a change occurs in a system. One type of such changes in climate
may be due to localized sources in space and time like a volcanic eruption.
The sulfur dioxide aerosols injected into the stratosphere by volcanic
eruptions can absorb incoming solar radiation, increasing temperatures in the
stratosphere and decreasing surface temperatures (Marshall et al.,, 2022). In
climate simulation experiments, volcanic eruptions lead to a variety of
changes in climate, including fewer temperature extremes (Paik & Min,, 2018)
and changes to the El Niño cycle (Khodri et al.,, 2017). Climate
geoengineering events, a broad class of potential human interventions on the
Earth’s climate, may also lead to changes in the climate (NASEM,, 2021; Lenton
& Vaughan,, 2009). Of these events, an injection of aerosols into the
stratosphere is the most similar to a volcanic eruption. Other interventions
include marine cloud brightening (increasing the reflectance of low-level
clouds) and removal of carbon dioxide from the atmosphere (Lenton & Vaughan,,
2009). The effects of such interventions are often not well understood, and
approaches for evaluating their downstream effects in the climate system are
needed when planning such interventions. Furthermore, methodology for
detecting such an intervention and evaluating its effect on the climate can
inform policy and security interests, since such interventions may
disproportionately affect different countries and could be done unilaterally
(NASEM,, 2021).

In this work, we propose and evaluate statistical methodology in the context
of the June 1991 eruption of Mt. Pinatubo in the Philippines. The eruption
injected approximately 20 million tons of sulfur dioxide in the stratosphere,
which spread over the tropics in the following weeks (Bluth et al.,, 1992).
The eruption has been linked using observations to increases in stratospheric
temperature (Labitzke & McCormick,, 1992) and decreases in surface
temperatures (Parker et al.,, 1996). The eruption size and the availability of
atmospheric observations at the time make it an attractive test case for
statistical methodology evaluating climate impacts from localized sources. We
are interested in a few questions. Can we detect changes in the climate
associated with Mt. Pinatubo in the presence of complex climate variability?
What is that heterogeneous effect on the climate? What does this test case
tell us about our ability to detect and evaluate climate geoengineering?

Various approaches have been proposed for detecting changes in climate data.
The work Reeves et al., (2007) compares approaches for basic univariate time
series. More recently, Lund et al., (2022) discusses important aspects and
challenges of changepoint detection in climate time series data, including
controlling for seasonality. Both reviews focus on time series and ignore the
potential spatial nature of the data. Majumdar et al., (2005) considers the
spatial problem yet estimates a change time that is shared across locations.
One work that addresses the spatial problem and heterogeneity across locations
is Moradi et al., (2023), though adjusting for seasonality is not addressed.

In the statistics literature, one perspective for changepoint estimation is
functional time series (for example, Hörmann & Kokoszka,, 2012), where we wish
to detect a distributional change in a sequence of functions. In Figure 1, we
present stratospheric temperature data from the MERRA-2 climate reanalysis at
one location as a functional time series (GMAO,, 2015).

Figure 1: MERRA-2 climate reanalysis data presented as a functional time
series (GMAO,, 2015) at $115^{\circ}$E and $10^{\circ}$N.

In this setting, we let $k=1,\dots,N$ denote the year and observe a function
of temperature $Y_{k}(u)$ of $u\in[0,1]$ where $u\cdot 365.25$ is the day of
the year. The functional time series at-most-one-change (AMOC) model is

$\displaystyle Y_{k}(u)$
$\displaystyle=\mu(u)+\mathbbm{1}(k>\tau)\delta(u)+\epsilon_{k}(u).$

Here $\mu(u)$ is the mean function of the data before the changepoint $\tau$,
after which the mean is $\mu(u)+\delta(u)$. The error $\epsilon_{k}(u)$ is
assumed to have zero mean and may be dependent in $u$ and $k$. Approaches for
detecting whether a change exists and for estimating $\tau$ have been proposed
in various works including Berkes et al., (2009); Aue et al., (2009, 2018);
Aston & Kirch, (2012), and we review these in more detail in Section 2.1.

In this work, we consider the problem where we observe such functional time
series at many spatial locations, and we aim to evaluate changes at each
location. That is, instead of looking at data at only one location as in
Figure 1, we are interested in heterogenous effects of Mt. Pinatubo at
different locations across the globe. There are some opportunities and
challenges with this additional structure to the problem. Spatial coherence in
the changepoints should allow one to borrow strength from different locations
to improve their estimates. There may also be spatial dependence in nuisance
variability which challenges our ability to borrow this strength. The full
heterogenous model will allow for a spatially-varying time and size of a
detected change. Inferring changepoints for this model is a challenging
multiple testing problem with null hypotheses $H_{0i}:\\{\textrm{no change at
location }\boldsymbol{s}_{i}\\}$ for $i=1$, $2$, $\dots$, $n$, where we
observe data at $n$ different locations indexed by $\boldsymbol{s}_{i}$.

Changepoint detection for spatially-dependent functional time series began in
Gromenko et al., (2017) which was later extended in Zhao et al., (2019) and
Dette & Quanz, (2022). While these approaches address spatial dependence in
the functional data, one changepoint time for all locations is assumed (as in
Majumdar et al.,, 2005), resulting in no spatial heterogeneity and a single
determination if a changepoint occurred. This may be appropriate for data
defined on small spatial regions, but it is reductive for globally-defined
data. We are particularly inspired by the model recently proposed by Wang et
al., (2022). It has many attractive components mentioned earlier: a spatially
varying changepoint time $\tau(\boldsymbol{s})$, spatially-varying mean
function $\mu(\boldsymbol{s},u)$ and change function
$\delta(\boldsymbol{s},u)$, and spatially-dependent error
$\epsilon_{k}(\boldsymbol{s},u)$. For the data $Y_{k}(\boldsymbol{s}_{i},u)$
now indexed by $\boldsymbol{s}$, the model is

$\displaystyle Y_{k}(\boldsymbol{s}_{i},u)$
$\displaystyle=\mu(\boldsymbol{s}_{i},u)+\mathbb{I}(k>\tau(\boldsymbol{s}_{i}))\delta(\boldsymbol{s}_{i},u)+\epsilon_{k}(\boldsymbol{s}_{i},u).$
(1)

We also assume that we observe the function $Y_{k}(\boldsymbol{s}_{i},u)$ at
points $u_{j}$, $j=1,\dots,m$, with $u_{j}$ evenly spaced on $[0,1]$. Wang et
al., (2022) propose a Bayesian hierarchical model that vastly improves
estimates of $\tau(\boldsymbol{s})$ compared to approaches that ignore space.

We will propose an approach for estimation of the model in (1) that is more
broad, applicable, and direct compared to Wang et al., (2022). For example,
Wang et al., (2022) consider a single family of test statistics called “fully-
functional,” and we will develop approaches for both “fully-functional” and
“score-based” statistics. While Wang et al., (2022) focus on the spatial
estimation of $\tau(\boldsymbol{s})$, we also show that one can improve
changepoint detection using spatial information. One major contrast is that,
while Wang et al., (2022) model the changepoint test statistics, we directly
model the original data. There are benefits and drawbacks for each approach.
For example, to facilitate modeling the original data, we assume its
normality, a much more restrictive condition than considered in Wang et al.,
(2022). Conversely, modeling the test statistics poses its own challenges.
Wang et al., (2022) use a piecewise-linear and normal approximation to
nonnegative test statistics; while they show that this works in practice, such
a formulation is somewhat unnatural. This also results in a relatively
complicated and computationally-intensive Bayesian hierarchical model. We will
show that by modeling the dependence in the original data, we can improve both
changepoint detection and estimates of the changepoint time using relatively
simple tools for estimating the model (1) without proposing a new model for
the test statistics.

We evaluate our approach in the context of a statistical simulation study, a
simplified climate simulation, and climate reanalysis, which represent
successively more challenging applications. The simulation study shows that,
by sharing information across space, one can dramatically increase power for
changepoint detection, use multiple-testing corrections to control the false
discovery rate, and estimate spatially-varying changepoints with more
accuracy. We next apply changepoint techniques to simplified climate model
output based on the Held-Suarez-Williamson model (HSW, Held & Suarez,, 1994;
Williamson et al.,, 1998) in Hollowed, (2022) (referred to as HSW++). This
simulates a Mt. Pinatubo-like injection of aerosols in the stratosphere with
simplified climate processes. With multiple-testing corrections on spatially-
predicted test statistics, we detect spatially-varying changes in
stratospheric temperature in line with expectations. The analysis on climate
reanalysis data (MERRA-2, GMAO,, 2015) is the most challenging setting in this
work. When using spatial information, we detect stratospheric temperature
increases and evaluate them in the context of Mt. Pinatubo by comparing
estimates of $\tau(\boldsymbol{s})$ and $\delta(\boldsymbol{s},u)$. Overall,
this work establishes the spatially-correlated functional time series
changepoint detection as a useful tool to detect changes in the climate due to
localized sources.

We discuss the organization of the rest of the paper. In Section 2, we
introduce the changepoint detection procedures used, including the proposed
spatially-predicted test statistics. In Section 3, we propose a basic
estimation procedure for a spatially-dependent functional data model indexed
on the sphere (written as $\mathbb{S}^{2}$), which to our knowledge has not
yet been addressed. Section 4 details the simulation study, and Section 5
presents the two data analyses. We discuss areas of improvement for future
analyses in Section 6. Accompanying the paper is a supplement with additional
details.

##  2 Changepoint detection for functional time series

###  2.1 Application at each location separately

Our baseline detection procedure computes test statistics using data only at
one location at a time, ignoring spatial coherency and dependence in the data.
We review this literature, focusing on “score-based” and “fully-functional”
statistics. Other approaches for detecting changes in functional data include
the fused-lasso approach of Harris et al., (2022).

####  2.1.1 Score-based approaches

The score-based group of detection procedures was proposed in Berkes et al.,
(2009) and Aue et al., (2009). This approach estimates a model under the null
hypothesis and uses functional principal component scores to evaluate a
change. That is, one estimates the model at location $\boldsymbol{s}_{i}$ of
$Y_{k}(\boldsymbol{s}_{i},u_{j})=\mu^{\boldsymbol{s}_{i}}(u_{j})+\sum_{q=1}^{Q}Z^{\boldsymbol{s}_{i}}_{qk}\phi_{q}^{\boldsymbol{s}_{i}}(u_{j})$
for some fixed $Q$. We refer to
$\\{\phi_{q}^{\boldsymbol{s}_{i}}(\cdot)\\}_{q}$ and
$\\{Z_{qk}^{\boldsymbol{s}_{i}}\\}_{q}$ as the functional principal components
and scores of the data at $\boldsymbol{s}_{i}$, respectively. Here, the
superscript $\boldsymbol{s}_{i}$ indicates that the estimates only use the
data at $\boldsymbol{s}_{i}$. To test the null hypothesis of
$\delta(\boldsymbol{s}_{i},u)=0$ for all $u\in[0,1]$, one may use the test
statistic

$\displaystyle
T_{score}^{\boldsymbol{s}_{i}}=\frac{1}{N^{2}}\sum_{q=1}^{Q}\frac{1}{\lambda_{q}^{\boldsymbol{s}_{i}}}\sum_{t=1}^{N}\left(\sum_{k=1}^{t}Z_{qk}^{\boldsymbol{s}_{i}}-\frac{t}{N}\sum_{k=1}^{N}Z_{qk}^{\boldsymbol{s}_{i}}\right)^{2},$
(2)

where
$\lambda_{q}^{\boldsymbol{s}_{i}}=\textrm{Var}(Z^{\boldsymbol{s}_{i}}_{qk})$
is the variance of the $q$-th principal component at location
$\boldsymbol{s}_{i}$. The test statistic is a normalized sum of CUSUM
(cumulative sum) statistics of the first $Q$ scores. Under conditions on the
error process $\epsilon_{k}(\boldsymbol{s}_{i},u)$, Berkes et al., (2009) show
that under the null hypothesis
$T_{score}^{\boldsymbol{s}_{i}}\overset{d}{\to}\sum_{q=1}^{Q}\int_{0}^{1}B_{q}^{2}(x)dx$
as $N\to\infty$ where $\\{B_{q}(\cdot)\\}_{q=1}^{Q}$ are independent, standard
Brownian bridges and $\overset{d}{\to}$ represents convergence in
distribution. Furthermore, they establish that if
$\delta(\boldsymbol{s}_{i},\cdot)$ is not orthogonal to the subspace spanned
by $\\{\phi_{q}^{\boldsymbol{s}_{i}}\\}_{q=1}^{Q}$, then
$T_{score}^{\boldsymbol{s}_{i}}$ converges in probability to $\infty$ as
$N\to\infty$. To estimate a changepoint time, Aue et al., (2009) propose and
provide asymptotic analysis of the estimator

$\displaystyle\hat{\tau}^{\boldsymbol{s}_{i}}_{score}$
$\displaystyle=\operatorname*{arg\,max}_{t=1,\dots,N}\frac{1}{N}\sum_{q=1}^{Q}\frac{1}{\lambda_{q}^{\boldsymbol{s}_{i}}}\left(\sum_{k=1}^{t}Z_{qk}^{\boldsymbol{s}_{i}}-\frac{t}{N}\sum_{k=1}^{N}Z_{qk}^{\boldsymbol{s}_{i}}\right)^{2}.$

Various work has expanded upon this score-based approach. For example, Jiao et
al., 2022a  aim to find a favorable set of basis functions for projection that
may not necessarily be eigenfunctions. Stoehr et al., (2021) and Jiao et al.,
2022b  propose detection procedures for changes in the functional covariance.
Chiou et al., (2019) use score-based test statistics to find multiple
changepoints. Aston & Kirch, (2012) analyze a “epidemic”-type changepoint,
where the mean later returns to its original level after a second time point.

####  2.1.2 Fully-functional approaches

A newer approach for changepoint detection for functional data was introduced
in Aue et al., (2018). This “fully-functional” test statistic relies upon the
raw data instead of estimated principal component scores, but the test
statistic takes a similar CUSUM-like form. Aue et al., (2018) introduce the
test statistic of $T_{ff}^{\boldsymbol{s}_{i}}=\max_{t}\lVert
S_{t,ff}^{\boldsymbol{s}_{i}}(\cdot)\rVert^{2}$, where

$\displaystyle
S_{t,ff}^{\boldsymbol{s}_{i}}(u)=\frac{1}{\sqrt{N}}\left\\{\sum_{k=1}^{t}Y_{k}^{\boldsymbol{s}_{i}}(u)-\frac{t}{N}\sum_{k=1}^{N}Y_{k}^{\boldsymbol{s}_{i}}(u)\right\\},\
t=0,\dots,N,$ (3)

$Y_{k}^{\boldsymbol{s}_{i}}(u)=Y_{k}(\boldsymbol{s}_{i},u)$ is viewed as data
at only that location, and $\lVert f\rVert^{2}=\int_{0}^{1}f^{2}(u)du$. To
estimate the distribution of $T_{ff}^{\boldsymbol{s}_{i}}$ under the null
hypothesis, one must estimate the long-run variance of the leading principal
components. These principal component variances
$\\{\lambda_{q,ff}^{\boldsymbol{s}_{i}}\\}$ are estimated assuming the
alternative hypothesis is true. In particular, it is shown that
$T_{ff}^{\boldsymbol{s}_{i}}\overset{d}{\to}\displaystyle\sup_{0\leq x\leq
1}\sum_{q=1}^{\infty}\lambda_{q,ff}^{\boldsymbol{s}_{i}}B_{q}^{2}(x)$ where
$\\{B_{q}(x)\\}$ again are independent standard Brownian bridges. In practice,
the infinite sum must be truncated and approximated with the first $Q$ terms.
To estimate $\tau(\boldsymbol{s}_{i})$, a natural estimator is
$\hat{\tau}_{ff}^{\boldsymbol{s}_{i}}=\operatorname*{arg\,max}_{t}\lVert
S_{t,ff}^{\boldsymbol{s}_{i}}(\cdot)\rVert^{2},$ and Aue et al., (2018) also
justify this estimator in theoretical results.

The fully-functional approach has also spawned further research. Horváth et
al., (2022), Aue et al., (2020), and Dette & Kutta, (2021) detect changes in
the covariance function using fully-functional or similar test statistics.
Rice & Zhang, (2022) estimate multiple changepoints, and Wegner & Wendler,
(2022) propose a robust fully-functional test statistic. Tucker & Yarger,
(2022) develop test statistics for functional data with phase variability.

###  2.2 Assumed model for error

To share information across locations, some structure in the error process
$\epsilon_{k}(\boldsymbol{s},u)$ is assumed, which influences how one analyzes
the data. Most critically, we will incorporate spatial dependence in
$\epsilon_{k}(\boldsymbol{s},u)$ through its principal component scores. We
assume that

$\displaystyle\epsilon_{k}(\boldsymbol{s}_{i},u_{j})$
$\displaystyle=\sum_{q=1}^{Q}Z_{qk}(\boldsymbol{s}_{i})\phi_{q}(u_{j})+W_{ijk},$

where $W_{ijk}\sim\mathcal{N}(0,\xi^{2})$ is normally-distributed measurement
error with variance $\xi^{2}>0$, $\\{\phi_{q}(\cdot)\\}$ are a set of
principal component functions, and $\\{Z_{qk}(\boldsymbol{s})\\}$ are spatial
fields of principal component scores. The scores are further decomposed
$Z_{qk}(\boldsymbol{s}_{i})=\tilde{Z}_{qk}(\boldsymbol{s}_{i})+U_{iqk}$,
where, for each $q$ and $k$, $\tilde{Z}_{qk}(\boldsymbol{s}_{i})$ is a
spatially-dependent random field with a Matérn covariance function, and
$\\{U_{iqk}\\}$ are independent across all $i$, $q$, and $k$ and distributed
$\mathcal{N}(0,\gamma_{q}^{2}(\boldsymbol{s}))$. As a simplified example, if
the smoothness of the Matérn covariance is set to $1/2$, we obtain

$\displaystyle\begin{split}\textrm{Cov}(Z_{q_{1}k_{1}}(\boldsymbol{s}_{1}),Z_{q_{2}k_{2}}(\boldsymbol{s}_{2}))&=\mathbbm{1}(q_{1}=q_{2},k_{1}=k_{2})\sigma_{q_{1}}(\boldsymbol{s}_{1})\sigma_{q_{1}}(\boldsymbol{s}_{2})\textrm{exp}\left(-\frac{\left\lVert\boldsymbol{s}_{1}-\boldsymbol{s}_{2}\right\rVert}{\alpha_{q_{1}}}\right)\\\
&~{}~{}~{}~{}+\mathbbm{1}(q_{1}=q_{2},k_{1}=k_{2},\boldsymbol{s}_{1}=\boldsymbol{s}_{2})\gamma_{q_{1}}^{2}(\boldsymbol{s}_{1}).\end{split}$
(4)

Here $\sigma_{q_{1}}^{2}(\boldsymbol{s}_{1})$ and
$\gamma_{q_{1}}^{2}(\boldsymbol{s}_{1})$ are spatially-varying variance terms
and $\alpha_{q}$ is a range or scale parameter describing the decay of
dependence. The spatially-varying variance term is based on Stein, (2005)
generalizing the approach of Paciorek & Schervish, (2006) that introduces
nonstationary forms for the parameter $\alpha_{q}$. We assume that
$\\{W_{ijk}\\}$, $\\{U_{iqk}\\}$, and
$\\{\tilde{Z}_{qk}(\boldsymbol{s}_{i})\\}$ are jointly normal and independent
from each other. A discussion and comparison of these assumptions is included
in Section S1 of the supplement. We also evaluate the normality assumption in
Section S5.

###  2.3 An approach for spatially-dependent functional time series

We next introduce our approach for improving upon test statistics using
spatial information, consisting of three steps detailed in the following three
subsections. The motivation for spatial prediction is to reduce the
independent error in the principal component scores $U_{iqk}$ or measurement
error $W_{ijk}$. We will focus on data observed on a fine grid, so spatial
interpolation is less of a motivation. Also, it is unclear if one should
estimate the model under the null or alternative hypothesis. We will evaluate
both in our simulation study.

####  2.3.1 Estimate spatio-functional model

We outline our approach for estimating the entire model (1) under the null or
an alternative hypothesis, resulting in the estimate
$\hat{\mu}(\boldsymbol{s},u)$, possibly $\hat{\delta}(\boldsymbol{s},u)$, and
an estimated covariance structure of $\epsilon_{k}(\boldsymbol{s},u)$. To
estimate a model under an alternative hypothesis, one should first specify a
“pilot” estimate of $\hat{\tau}(\boldsymbol{s}_{i})$. For data at one
location, one can estimate each of the $N$ models taking
$\hat{\tau}(\boldsymbol{s}_{1})\in\\{1,2,\dots,N\\}$. With spatial dependence,
there are $N^{n}$ models, and one cannot estimate all of them. Sensible pilot
estimates of $[\hat{\tau}(\boldsymbol{s}_{i})]_{i=1}^{n}$ are those in Section
2.1. Assuming the null hypothesis requires no pilot estimate as
$\hat{\tau}(\boldsymbol{s}_{i})=N$ for all $i$.

First, the fixed effects $\mu(\boldsymbol{s},u)$ and
$\delta(\boldsymbol{s},u)$ are estimated. We use penalized least squares
estimators for these terms. Exact formulation of estimators may depend on the
domain of $\boldsymbol{s}$ and $u$. In Section 3.1, we specifically outline an
approach where $\boldsymbol{s}\in\mathbb{S}^{2}$ and $u\in[0,1]$.

We use a similar approach to estimate the marginal covariance structure in
$u$. If we define
$\boldsymbol{Y}_{k}(\boldsymbol{s}_{i})=\left[Y_{k}(\boldsymbol{s}_{i},u_{j})\right]_{j=1}^{m}$
and
$\hat{\boldsymbol{\mu}}_{k}(\boldsymbol{s}_{i})=\left[\hat{\mu}(\boldsymbol{s}_{i},u_{j})+\mathbbm{1}(k>\hat{\tau}(\boldsymbol{s}_{i}))\hat{\delta}(\boldsymbol{s}_{i},u_{j})\right]_{j=1}^{m}$,
an empirical covariance estimate in $u$ collapsed over all years and spatial
locations is

$\displaystyle\hat{\boldsymbol{C}}_{u}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{N}\sum_{k=1}^{N}\left(\boldsymbol{Y}_{k}(\boldsymbol{s}_{i})-\hat{\boldsymbol{\mu}}_{k}(\boldsymbol{s}_{i})\right)\left(\boldsymbol{Y}_{k}(\boldsymbol{s}_{i})-\hat{\boldsymbol{\mu}}_{k}(\boldsymbol{s}_{i})\right)^{\top}.$

In practice we include a smoothing penalty on the resulting covariance
estimate, and we implement generalized cross-validation (c.f. Wahba,, 1990) to
choose smoothing parameters.

Let $\hat{\boldsymbol{\Phi}}\in\mathbb{R}^{m\times Q}$ be the matrix of the
first $Q$ functional principal components based on $\hat{\boldsymbol{C}}_{u}$
with $Q<m$ evaluated at $[u_{j}]_{j=1}^{m}$. The principal component scores
are estimated with
$[Z_{qk}(\boldsymbol{s}_{i})]_{q=1}^{Q}=\hat{\boldsymbol{\Phi}}^{\top}\left(\boldsymbol{Y}_{k}(\boldsymbol{s}_{i})-\hat{\boldsymbol{\mu}}_{k}(\boldsymbol{s}_{i})\right)$
omitting the “hat” on $Z_{qk}(\boldsymbol{s})$ for notational flexibility
later. We then form estimators $\hat{\gamma}_{q}^{2}(\boldsymbol{s})$ and
$\hat{\sigma}_{q}^{2}(\boldsymbol{s})$ based on Remark S2 of the supplement.

We next turn to the estimation of the dependence structure of
$\\{Z_{qk}(\boldsymbol{s})\\}$, for example, the parameters
$\\{\alpha_{q}\\}_{q=1}^{Q}$ and smoothness parameters
$\\{\nu_{q}\\}_{q=1}^{Q}$ assumed constant over $k$. Let
$\boldsymbol{\theta}_{q}=(\alpha_{q},\nu_{q},\hat{\gamma}_{q}^{2}(\boldsymbol{s}),\hat{\sigma}_{q}^{2}(\boldsymbol{s}))$
denote the covariance’s parameters, yet we treat
$\hat{\gamma}_{q}^{2}(\boldsymbol{s})$ and
$\hat{\sigma}_{q}^{2}(\boldsymbol{s})$ as fixed since their flexible structure
is estimated previously. Let
$\boldsymbol{Z}_{qk}=[Z_{qk}(\boldsymbol{s}_{i})]_{i=1}^{n}$ and
$\boldsymbol{C}_{s}(\boldsymbol{\theta}_{q})=\textrm{Cov}(\boldsymbol{Z}_{qk},\boldsymbol{Z}_{qk}|\boldsymbol{\theta}_{q})$
be its covariance matrix using $\boldsymbol{\theta}_{q}$ and specified by a
form like (4). Then, the Gaussian likelihood for the $q$-th principal
component is

$\displaystyle\ell_{q}(\boldsymbol{\theta}_{q})$
$\displaystyle=-\frac{1}{2}\sum_{k=1}^{N}\left(\log(\textrm{det}(\boldsymbol{C}_{s}(\boldsymbol{\theta}_{q})))+\boldsymbol{Z}_{qk}^{\top}\boldsymbol{C}_{s}(\boldsymbol{\theta}_{q})^{-1}\boldsymbol{Z}_{qk}+n\log(2\pi)\right).$

In our simulation study, we directly evaluate the likelihood, yet for the data
application some approximation is necessary as detailed in Section 3. Maximum
likelihood estimation through numerical optimization can then proceed
separately for each $q=1,\dots,Q$.

####  2.3.2 Use model for prediction

Referring to the notation in Section 2.2, our goal is to form an estimate of
$\tilde{\boldsymbol{Z}}_{qk}=\left[\tilde{Z}_{qk}(\boldsymbol{s}_{i})\right]_{i=1}^{n}$,
a spatially-smoothed version of the scores $\boldsymbol{Z}_{qk}$. Based on the
multivariate normal distribution, the conditional expectation using estimated
parameters $\hat{\boldsymbol{\theta}}_{q}$ for each $q$ and $k$ is

$\displaystyle\hat{\boldsymbol{Z}}_{qk}=\left[\hat{Z}_{qk}(\boldsymbol{s}_{i})\right]_{i=1}^{n}=\mathbb{E}\left[\tilde{\boldsymbol{Z}}_{qk}\middle|\boldsymbol{Z}_{qk},\hat{\boldsymbol{\theta}}_{q}\right]$
$\displaystyle=\left[\tilde{\boldsymbol{C}}_{s}\left(\hat{\boldsymbol{\theta}}_{q}\right)\right]\left[\boldsymbol{C}_{s}\left(\hat{\boldsymbol{\theta}}_{q}\right)\right]^{-1}\boldsymbol{Z}_{qk}.$

Here, the $n\times n$ matrix
$\tilde{\boldsymbol{C}}_{s}\left(\hat{\boldsymbol{\theta}}_{q}\right)=\textrm{Cov}\left(\tilde{\boldsymbol{Z}}_{qk},\boldsymbol{Z}_{qk}\middle|\hat{\boldsymbol{\theta}}_{q}\right)$
uses the estimated covariance function. Then, one can form predictions of the
original data with reduced variance

$\displaystyle\hat{Y}_{k}(\boldsymbol{s}_{i},u_{j})$
$\displaystyle=\hat{\mu}(\boldsymbol{s}_{i},u_{j})+\mathbbm{1}(k>\hat{\tau}(\boldsymbol{s}_{i}))\hat{\delta}(\boldsymbol{s}_{i},u_{j})+\sum_{q=1}^{Q}\hat{Z}_{qk}(\boldsymbol{s}_{i})\hat{\phi}_{q}(u_{j}).$

####  2.3.3 Recompute test statistics using predictions

Using the predicted scores $\left\\{\hat{Z}_{qk}(\boldsymbol{s})\right\\}$ and
predicted data $\left\\{\hat{Y}_{k}(\boldsymbol{s}_{i},u_{j})\right\\}$, we
now recompute test statistics. Our primary score-based predicted test
statistic at location $\boldsymbol{s}_{i}$ is

$\displaystyle\hat{T}_{score}(\boldsymbol{s}_{i})=\frac{1}{N^{2}}\sum_{q=1}^{Q}\frac{1}{\hat{\lambda}^{*}_{q}(\boldsymbol{s}_{i})}\sum_{t=1}^{N}\left(\sum_{k=1}^{t}\hat{Z}_{qk}(\boldsymbol{s}_{i})-\frac{t}{N}\sum_{k=1}^{N}\hat{Z}_{qk}(\boldsymbol{s}_{i})\right)^{2}.$
(5)

where
$\hat{\lambda}^{*}_{q}(\boldsymbol{s}_{i})=\textrm{Var}\left(\hat{Z}_{qk}(\boldsymbol{s}_{i})\right)$
is a local-regression estimate of the marginal variance of the predicted
scores. We consider two additional flavors of a predicted score-based test
statistic. The first replaces $\hat{\lambda}^{*}_{q}(\boldsymbol{s}_{i})$ with
the estimate
$\hat{\lambda}_{q}(\boldsymbol{s}_{i})=\textrm{Var}(Z_{qk}(\boldsymbol{s}_{i}))$
defined in Remark S2; we say this has “unadjusted variances” since it uses the
variance of the original scores rather than the predicted scores. In addition,
one can recompute score-based test statistics at each location using the
predicted data $\left\\{\hat{Y}_{k}(\boldsymbol{s}_{i},u_{j})\right\\}$; we
refer to this approach as “recomputing individual test statistics.” For score-
based test statistics, we use predictions of scores and data assuming the null
hypothesis is true; the predicted scores can still be appropriately used
directly in (5) to detect a change. Otherwise, the change may be partitioned
in
$\hat{\mu}(\boldsymbol{s}_{i},u_{j})+\mathbbm{1}(k>\hat{\tau}(\boldsymbol{s}_{i}))\hat{\delta}(\boldsymbol{s}_{i},u_{j})$.
For fully functional test statistics, we use the predictions
$\left\\{\hat{Y}_{k}(\boldsymbol{s}_{i},u_{j})\right\\}$ to form and evaluate
test statistics.

###  2.4 Multiple testing

Inferring changes at multiple locations is fundamentally a multiple testing
problem. Throughout our simulation and data analysis, we will evaluate both
unadjusted for multiple testing and unadjusted for multiple testing p-values.
We adjust p-values using the Bonferroni and Benjamini-Hochberg procedures to
control family-wise error rate and false discovery rate, respectively
(Benjamini & Hochberg,, 1995). These procedures do not inherently consider the
spatial nature of the data, yet using them on predicted data may improve
results compared to using them without prediction. There is also substantial
literature on spatial multiple testing procedures, in the setting of
predetermined spatial clusters (Benjamini & Heller,, 2007), a Bayesian
framework (Risser et al.,, 2019), and an EM framework (Yun et al.,, 2022).
Since multiple testing is not the main focus of this work, we evaluate the
more simple Bonferroni and Benjamini-Hochberg approaches.

##  3 Spherically-indexed functional data framework

Here, we introduce our framework to estimate the functional data model (1)
when the spatial information $\boldsymbol{s}$ is indexed on the unit sphere
$\mathbb{S}^{2}$, which is used to approximate the surface of the Earth in our
analysis. In the spatial statistics community, there has been substantial
interest in modeling random fields on $\mathbb{S}^{2}$. See Jeong et al.,
(2017) for a review, and more recent contributions include Porcu et al.,
(2019) and Alegría et al., (2020). Modeling random fields on $\mathbb{S}^{2}$
can be more challenging in some aspects compared to modeling those on
$\mathbb{R}^{d}$. For example, the choice of distance (i.e. chordal or
geodesic distance) on $\mathbb{S}^{2}$ is thoroughly debated (e.g. Krivoruchko
& Gribov,, 2020). Axial symmetric models (Jeong et al.,, 2017) are an example
of models specific to $\mathbb{S}^{2}$.

We are interested in modeling random fields of functional data
$Y_{k}(\boldsymbol{s},u)$ indexed on the sphere with
$\boldsymbol{s}\in\mathbb{S}^{2}$ and $u\in[0,1]$. Though we are unaware of
previous literature on this, many have extended random fields on a sphere to
more complicated settings. For example, modeling sphere-time processes has
received attention, though not from the functional-data perspective; see Porcu
et al., (2018). Other work treats the sphere $\mathbb{S}^{2}$ as the
functional domain of the data (see Caponera et al.,, 2022, and references
within) or treats trajectories on $\mathbb{S}^{2}$ as functional data (for
example, Dai & Müller,, 2018; Su et al.,, 2014).

As we build our approach, we have a few guiding principles. For modeling the
principal component scores, a nonstationary model is necessary, especially
with respect to their variances. We also aim for an implementation that is
possible and straightforward on the scale of global data with potentially tens
of thousands of locations. We base the model for each principal component
score on Cao et al., (2022) to reach these goals. We first discuss the
estimation of the mean and change functions.

###  3.1 Estimation of fixed effects

To estimate $\mu(\boldsymbol{s},u)$ and $\delta(\boldsymbol{s},u)$, we propose
a standard penalized optimization-based estimator. For illustration, we focus
on $\mu(\boldsymbol{s},u)$, and estimate

$\displaystyle\hat{\mu}(\boldsymbol{s},u)=\operatorname*{arg\,min}_{\beta}\sum_{i=1}^{n}\sum_{j=1}^{m}\left(\overline{Y}^{\mu}(\boldsymbol{s}_{i},u_{j})-\beta(\boldsymbol{s}_{i},u_{j})\right)^{2}+\textrm{Pen}(\zeta).$

where
$\overline{Y}^{\mu}(\boldsymbol{s}_{i},u)=\textrm{mean}\left\\{Y_{k}(\boldsymbol{s}_{i},u)\textrm{
such that }k\leq\hat{\tau}(\boldsymbol{s})\right\\}$ is the empirical average
estimate based on an estimate $\hat{\tau}(\boldsymbol{s})$, and
$\textrm{Pen}(\zeta)$ is a smoothness penalty on the solution $\beta$ with
smoothing parameter $\zeta>0$. A similar estimator is used to estimate
$\delta(\boldsymbol{s},u)$. We represent $\beta(\boldsymbol{s},u)$ using a
product basis of spherical harmonic functions (of $\boldsymbol{s}$, sphunif,
García-Portugués & Verdebout,, 2021, 2018) and cubic B-splines (of $u$). For
the penalty $\textrm{Pen}(\zeta)$, we apply smoothing with respect to $u$
only. Since spherical harmonic functions are not local, using a large basis
can be computationally infeasible, so we use instead a regression approach
similar to Jeong et al., (2017). We provide details in Section S2.

###  3.2 Estimation of spatio-functional dependence structure

We focus on modeling the spatial dependence in principal component scores
$Z_{qk}(\boldsymbol{s}_{i})$; mostly following the notation of Cao et al.,
(2022), we take

$\displaystyle\textrm{Cov}(Z_{qk}(\boldsymbol{s}_{1}),Z_{qk}(\boldsymbol{s}_{2})|\boldsymbol{\theta}_{q})$
$\displaystyle=\sigma_{q}(\boldsymbol{s}_{1})\sigma_{q}(\boldsymbol{s}_{2})c_{q}(\boldsymbol{s}_{1},\boldsymbol{s}_{2})\rho_{q}(d_{q}(\boldsymbol{s}_{1},\boldsymbol{s}_{2}))+\mathbbm{1}(\boldsymbol{s}_{1}=\boldsymbol{s}_{2})\gamma_{q}^{2}(\boldsymbol{s}_{1}),$

where $c_{q}(\boldsymbol{s}_{1},\boldsymbol{s}_{2})$ is a normalizing
constant, $\rho_{q}(\cdot)$ is an isotropic covariance function, and
$d_{q}(\boldsymbol{s}_{1},\boldsymbol{s}_{2})$ is a distance function. Here,
we take

$\displaystyle d_{q}(\boldsymbol{s}_{1},\boldsymbol{s}_{2})^{2}$
$\displaystyle=2(\boldsymbol{s}_{1}-\boldsymbol{s}_{2})^{\top}(\boldsymbol{\Sigma}_{q}(\boldsymbol{s}_{1})+\boldsymbol{\Sigma}_{q}(\boldsymbol{s}_{2}))^{-1}(\boldsymbol{s}_{1}-\boldsymbol{s}_{2})$
$\displaystyle c_{q}(\boldsymbol{s}_{1},\boldsymbol{s}_{2})^{2}$
$\displaystyle=|\boldsymbol{\Sigma}_{q}(\boldsymbol{s}_{1})|^{1/2}|\boldsymbol{\Sigma}_{q}(\boldsymbol{s}_{2})|^{1/2}|(\boldsymbol{\Sigma}_{q}(\boldsymbol{s}_{1})+\boldsymbol{\Sigma}_{q}(\boldsymbol{s}_{2}))/2|^{-1}.$

The matrix $\boldsymbol{\Sigma}_{q}(\boldsymbol{s})$ specifies the spatially-
varying anisotropy and depends on the angle parameter
$\kappa_{q}(\boldsymbol{s})$ and length scale parameters
$\chi_{q1}(\boldsymbol{s})$ and $\chi_{q2}(\boldsymbol{s})$. For the exact
construction of $\boldsymbol{\Sigma}_{q}(\boldsymbol{s})$ from these
parameters, we refer the reader to Cao et al., (2022). Following their
implementation, we take $\kappa_{q}(\boldsymbol{s})=\kappa_{q}$. For the
length scale parameters, we take
$\chi_{qr}(\boldsymbol{s})=\textrm{exp}\left(\beta_{qr0}+\beta_{qr1}\sin(l)+\beta_{qr2}L+\beta_{qr3}L^{2}\right)$,
for $r=1,2$ and $l$ is the longitude and $L$ is the latitude of
$\boldsymbol{s}$. This differs slightly from Cao et al., (2022) as we
introduce a quadratic term of latitude, allowing for behavior in parameters
symmetric about the equator. For each $\rho_{q}(\cdot)$, a Matérn covariance
function is used, each with a smoothness parameter $\nu_{q}$ also estimated.
Therefore, $\boldsymbol{\theta}_{q}=\\{\nu_{q}$, $\kappa_{q}$,
$\\{\beta_{qr\ell}\\}_{r=1,\ell=0}^{2,3}$, $\hat{\sigma}_{q}^{2}(\cdot)$,
$\hat{\gamma}_{q}^{2}(\cdot)\\}$.

This covariance is attractive because it allows for a spatially-varying
variance and decay parameters as well as control over the smoothness of the
process. This covariance is also implemented with a Vecchia approximation in
Cao et al., (2022), which makes implementation straightforward for analyses on
the scope considered in this work. Instead of working with and inverting the
entire dense covariance matrix
$\boldsymbol{C}_{s}(\boldsymbol{\theta}_{q})\in\mathbb{R}^{n\times n}$ for
each likelihood evaluation or prediction, this instead builds a sparse
representation of the precision matrix
$(\boldsymbol{C}_{s}(\boldsymbol{\theta}_{q}))^{-1}$ based on a neighborhood
structure (Katzfuss et al.,, 2020). We also avoid recomputations over $k$ when
working with $\\{Z_{qk}(\boldsymbol{s})\\}$ for all $k=1,\dots,N$. The
parameters $\boldsymbol{\theta}_{q}$ (except $\hat{\sigma}_{q}^{2}(\cdot)$ and
$\hat{\gamma}_{q}^{2}(\cdot)$) are estimated by maximum likelihood as detailed
in Section 2.3.1.

##  4 Simulation

In this section, we consider a simulation study which finds that spatially-
predicted versions of test statistics improve both the detection of
changepoints and the estimation of $\tau(\boldsymbol{s})$. For illustrative
use, we consider a spatial domain on a grid in $\mathbb{R}^{2}$.

###  4.1 Simulation design

We again take the model defined in (1). We consider $N=15$, $n=300$ locations
on a $15\times 20$ grid on $[0,1]\times[0,1]$, and $m=40$ measurements in $u$
per location evenly-spaced on $[0,1]$. To specify the fixed components, we
take, with $\boldsymbol{s}=(s_{1},s_{2})^{\top}$,

$\displaystyle\mu(\boldsymbol{s},u)$ $\displaystyle=\cos(\pi
u)\textrm{exp}\left(\frac{2(s_{1}-s_{2})}{2u+1}\right),$
$\displaystyle\delta(\boldsymbol{s},u)$
$\displaystyle=\underbrace{\eta}_{\textrm{signal
strength}}\underbrace{\left(us_{1}+u^{2}s_{2}-u^{3}(s_{1}+s_{2})\right)}_{\textrm{varying
change
function}}\underbrace{\mathbbm{1}(s_{1}-s_{2}>-0.3)(s_{1}-s_{2}+0.3)}_{\textrm{taper
near boundary}},$ $\displaystyle\tau(\boldsymbol{s})$
$\displaystyle=\frac{1}{4}\left\lceil
15+3\textrm{exp}(s_{1}+s_{2})\right\rceil\textrm{ if }s_{1}-s_{2}>-.3,\textrm{
and }N\textrm{ (no changepoint) otherwise},$

where $\lceil x\rceil$ is the ceiling function of $x$. We plot these functions
in Figure 2. The design of $\tau(\boldsymbol{s})$ and
$\delta(\boldsymbol{s},u)$ avoids sudden changes in space between locations
with or without a change, which is realistic in our application, since
circulation prevents hard spatial boundaries in stratospheric temperature.
Thus, the locations with a change near the boundary with “no change” are
hardest to detect. We vary the signal strength with
$\eta\in\\{0.5,2,6,10,14,18,22\\}$.

Figure 2: (Far Left) Original means and (Left) means after the changepoint for
selected locations, with $\eta=10$. (Right) Changepoint surface
$\tau(\boldsymbol{s})$, with gray tiles having no change. (Far Right) Change
function for a fixed $u=0.5$.

Turning to the specification of $\epsilon_{k}(\boldsymbol{s},u)$, we take
$Q=3$ principal components based on

$\displaystyle\epsilon_{k}(\boldsymbol{s}_{i},u_{j})$
$\displaystyle=Z_{1k}(\boldsymbol{s}_{i})+Z_{2k}(\boldsymbol{s}_{i})2^{-\frac{1}{2}}\sin(2\pi
u_{j})+Z_{3k}(\boldsymbol{s}_{i})2^{-\frac{1}{2}}\cos(2\pi u_{j})+W_{ijk}.$

We decompose
$Z_{qk}(\boldsymbol{s}_{i})=\tilde{Z}_{qk}(\boldsymbol{s}_{i})+U_{iqk}$, and
the variances of $\tilde{Z}_{qk}(\boldsymbol{s}_{i})$ are spatially-varying:
$\sigma_{1}(\boldsymbol{s})=|s_{1}|+|s_{2}|+0.01$,
$\sigma_{2}(\boldsymbol{s})=(25)^{-1}(3|s_{1}|+27|s_{2}|)+0.01$, and
$\sigma_{3}(\boldsymbol{s})=(25)^{-1}(9|s_{1}|+|s_{2}|)+0.01$, while
$\textrm{Var}(U_{iqk})=\gamma_{q}^{2}(\boldsymbol{s}_{i})=0.3$ for all $q$ and
$\boldsymbol{s}_{i}$. We take $\textrm{Var}(W_{ijk})=0.16$.

We consider two settings of the absence or presence of spatial dependence in
$\epsilon_{k}(\boldsymbol{s},u)$:

  * •

Independent model: We take
$\textrm{Cov}(Z_{q_{1}k_{1}}(\boldsymbol{s}_{1}),Z_{q_{2}k_{2}}(\boldsymbol{s}_{2}))=\left(\sigma_{q_{1}}^{2}(\boldsymbol{s}_{1})+\gamma_{q_{1}}^{2}(\boldsymbol{s}_{1})\right)$
if $q_{1}=q_{2}$, $k_{1}=k_{2}$, and $\boldsymbol{s}_{1}=\boldsymbol{s}_{2}$,
and $0$ otherwise. Here, $\sigma_{q}^{2}(\boldsymbol{s})$ and
$\gamma_{q}^{2}(\boldsymbol{s})$ are not individually identifiable. Using
Remark S2, $\hat{\sigma}_{q}^{2}(\boldsymbol{s})$ estimates $0$ while
$\hat{\gamma}_{q}^{2}(\boldsymbol{s})$ estimates
$\sigma_{q}^{2}(\boldsymbol{s})+\gamma_{q}^{2}(\boldsymbol{s})$.

  * •

Dependent model: We take the model defined in (4) with
$(\alpha_{1},\alpha_{2},\alpha_{3})=(0.4,0.3,0.6)$.

An example of simulated $Z_{qk}(\boldsymbol{s}_{i})$ under the dependent model
are shown in the left panel of Figure 3.

Figure 3: (Left) Simulated principal component scores $Z_{qk}(\boldsymbol{s})$
for $q=1$, $2$, and $3$ and $k=1$. (Middle Right) Plots of the mean function
before and after the changepoint at one location $\boldsymbol{s}=(1,0.643)$
and (Far Right) simulated data for that location with $\eta=10$ and
$\tau(\boldsymbol{s})=8$.

We plot an example of simulated data at one location in Figure 3. We complete
100 replicates for each signal strength and for the independent and dependent
model.

###  4.2 Approaches compared

We first consider a variety of score-based test statistics. Our baseline
approach for comparison is the score-based test statistic computed at each
location while ignoring data at all other locations as in (2). For our
predicted test statistics, we use predictions assuming the null hypothesis is
true. We also vary the number of estimated principal components, taking $Q=3$
(the “true” number), $Q=4$, and $Q=5$; we also present results for $Q=4$ for
the “unadjusted variances” and “recomputed” variations of the predicted test
statistics described in Section 2.3.3. When estimating the model under the
null hypothesis, the change function may be effectively absorbed into the
covariance structure, which motivates using $Q>3$ to ensure we capture both
the change and residual variability.

For fully-functional test statistics, our baseline is again the test statistic
computed at each location. We also consider predictions under the global null
and four alternative models using results from the individual fully-functional
test statistic at each location: (“Alt”) a p-value $<0.05$, (“BH”) a
Benjamini-Hochberg-corrected p-value $<0.05$, (“Bonf”) a Bonferroni-corrected
p-value $<0.05$, (“All”) for all locations regardless of p-value. These models
lie on a spectrum between “the null is true everywhere” and “the alternative
is true everywhere.” We also compute the test statistics of Gromenko et al.,
(2017) for each simulation. We also tried applying Wang et al., (2022), but we
were not able to obtain convergent MCMC results that improved the estimate of
$\tau(\boldsymbol{s})$ in this setting.

###  4.3 Changepoint testing performance

In our simulation, we evaluate detection performance, with full results
presented in Section S3. We find that spatially predicting score-based test
statistics dramatically increases the detection power with slight degradations
in Type I error. Also, applying multiple testing criteria leads to control of
false discovery rate. After correction, only spatially-predicted test
statistics have substantial power. Finally, we compare results with Gromenko
et al., (2017), which, by testing for a single changepoint, is an
oversimplification in this setting.

###  4.4 Changepoint estimation performance

In addition to improving the test performance of the estimators, we can also
dramatically improve the estimation of the time of change. We first examine
performance across dependence structures and simulation strength in Figure 4
in terms of root-mean-squared-error between the true and estimated
changepoint. Here, we compare results for all locations with a true
changepoint, regardless of whether a changepoint was detected, so that we can
evaluate this performance straightforwardly and separately from testing.

Figure 4: Estimation year performance across simulation setups (dependence
structure and signal strength). (Left) Score-based test statistics. (Right)
Fully-functional test statistics. The scale of the plots are the same to ease
comparison.

Generally, the spatially-predicted test statistics improve upon test
statistics computed individually at each location. While there is no
improvement for lower signal strengths, at mid-to-high signal strengths, the
predicted test statistics can cut RMSE in half in the spatially-independent
setting. In the spatially-dependent setting, there are more modest but
substantial improvements. As signal strength increases from high to higher,
the gap between spatial and individual test statistics decreases, as the
change time is more easily estimated using data at each location.

Comparing results across the independent and dependent settings, the test
statistics computed at individual locations perform the same with or without
spatial dependence. This result is expected since the marginal distribution at
each location is the same, with spatial dependence only affecting how
different locations relate to each other. For the predicted test statistics,
adding spatial dependence makes the problem fundamentally harder by adding
spatially-correlated noise.

We also compare the performance of the score-based and fully-functional test
statistics. When computed individually at each location, the score-based test
statistics with 4 or 5 principal components perform better than the fully-
functional statistics. This result mostly holds for spatially-predicted test
statistics also. We see that predictions under alternative hypotheses do not
improve the test statistics much compared to other predicted test statistics.
Only using three principal components is substantially worse than using four
or five, since the change may be filtered into the covariance when estimating
the model under the null. An example of the changepoint estimation is given in
Section S3.

##  5 Data analysis

###  5.1 Application to HSW++ data

For an initial test case for the methodology in the context of global climate
data, we apply and evaluate the methodology on a simplified atmospheric model
branded as HSW++. This model builds upon the framework of Held & Suarez,
(1994) and Williamson et al., (1998) in the thesis work of Hollowed, (2022).
The primary goal is to evaluate the effects of a stratospheric aerosol
injection of $\textrm{SO}_{\textrm{2}}$, sulfate, and ash, designed to
simulate a volcanic eruption. While HSW++ is considerably simplified compared
to Earth Systems Models (ESMs), it provides a middle ground between a
statistical simulation study and an application to ESM output or climate
reanalysis data. HSW++ assumes that land with no elevation covers the entire
earth, there is no precipitation or water vapor, and there are no diurnal or
seasonal cycles. See Hollowed, (2022) for a detailed description of the model.

Figure 5: (Left) Longitudinal and ensemble average of 5 HSW++ ensemble members
with stratospheric injection. (Right) Example of HSW++ data for each ensemble
member for $38^{\circ}\textrm{E}$ and $18^{\circ}\textrm{N}$. Unless otherwise
noted, the vertical line represents the time of injection.

We apply the developed methodology separately to a 5-member ensemble run of
HSW++ where each member had a stratospheric aerosol injection event and
different initial conditions, as well as a “counterfactual” member with no
stratospheric aerosol injection event. Each ensemble member consists of data
defined on a $2^{\circ}\times 2^{\circ}$ latitude-longitude grid (16,380
locations total) over the course of 1,200 days. The stratospheric injection
occurs on the 180th day. We use 6-hourly averaged data resulting in 4,800
total data points per location and ensemble. Figure 5 shows a summary and
example of the stratospheric (near-50mb) temperature data. From the lack of
cyclical effects in the data, we can flexibly define $m$ and $N$ with $m=120$
and $N=40$, so that the “years” $k$ can now be naturally interpreted as a
“month” of $30$ days. Throughout, we use $Q=10$ principal components,
spherical harmonics up to degree $10$ resulting in $121$ spatial basis
functions, a bandwidth of $h=400$ kilometers, and $8$ neighbors for the
Vecchia approximation. We use spatially-predicted test statistics under the
null hypothesis and a Benjamini-Hochberg correction to the p-values.

At first, we focused on the “at-most-one-change” (AMOC) model described by
(1). The results for spatially-predicted test statistics are presented on the
left of Figure 6. While the third ensemble member had many changepoints
detected around the time of injection (month 6), for other ensemble members
many changes are detected near the middle of the time period (month 18).
Comparing with Figure 5, it appears that the AMOC changepoint model picks up a
decrease in temperature associated with a “return-to-normal” rather than the
immediate impact of the injection. For this reason, we began to consider an
“epidemic”-type changepoint model studied in Aston & Kirch, (2012). This model
assumes replaces the indicator $\mathbb{I}(k>\tau(\boldsymbol{s}_{i}))$ in (1)
with
$\mathbb{I}(\tau_{1}(\boldsymbol{s}_{i})<k\leq\tau_{2}(\boldsymbol{s}_{i}))$.
Ideally, the change $\tau_{1}(\boldsymbol{s})$ represents the injection’s
initial impact, while $\tau_{2}(\boldsymbol{s})$ represents the “return-to-
normal” previously detected. We use score-based statistics from Aston & Kirch,
(2012):

$\displaystyle T_{score,epi}^{\boldsymbol{s}_{i}}$
$\displaystyle=\frac{1}{N^{3}}\sum_{q=1}^{Q}\frac{1}{\lambda_{q}^{\boldsymbol{s}_{i}}}\sum_{t_{1}=1}^{N}\sum_{t_{2}>t_{1}}^{N}\left(\sum_{k=t_{1}}^{t_{2}}Z_{qk}^{\boldsymbol{s}_{i}}-\frac{(t_{2}-t_{1}+1)}{N}\sum_{k=1}^{N}Z_{qk}^{\boldsymbol{s}_{i}}\right)^{2},$

with convergence under the null hypothesis to $\displaystyle
T_{score,epi}^{\boldsymbol{s}_{i}}\overset{d}{\to}\sum_{q=1}^{Q}\int_{0}^{1}\int_{0}^{y}(B_{q}(x)-B_{q}(y))^{2}dx\
dy.$ See Aston & Kirch, (2012) for more information. Our approach for
including spatial information is directly applied by using predicted versions
of principal component scores.

Figure 6: Comparison of AMOC (Left) and epidemic (Middle and Right)
changepoint model for HSW++. (Middle) For estimates of
$\tau_{1}(\boldsymbol{s})$; (Right) for estimates of
$\tau_{2}(\boldsymbol{s})$. For each latitude, we plot the proportion of
locations (of 180 different longitudes) that had a change detected during that
month using Benjamini-Hochberg-corrected p-values.

The results for the epidemic-type changepoint are plotted on the right of
Figure 6. Estimates of $\tau_{1}(\boldsymbol{s})$ are concentrated around the
stratospheric aerosol injection event, and estimates of
$\tau_{2}(\boldsymbol{s})$ correspond more closely to the changepoints
detected with the AMOC model. There is variability in the results between
ensembles. For ensemble 2, very few changepoints were detected after
correction, and for ensemble 3, the estimates for $\tau_{2}(\boldsymbol{s})$
are later compared to other ensembles. For the counterfactual ensemble member,
after Benjamini-Hochberg correction only few changepoints were detected using
the AMOC model, while no changepoints were detected for the epidemic model.
This suggests that the changepoint detection procedures work well under a
“null hypothesis” in this climate model.

Figure 7: HSW++ results. (Left) Number of ensembles with detected epidemic
change after (Top) and before (Bottom) spatial prediction. (Top Right) Number
of ensembles with detected epidemic change after Benjamini-Hochberg correction
on spatially-predicted test statistics. (Bottom Right) Estimated average
magnitude of change surface.

In Figure 7, we plot results on the latitude-longitude grid. On the left, we
compare the number of ensembles (of 5 with injection) that the locations had a
detected epidemic changepoint using individual test statistics and test
statistics based on spatial prediction, both without p-value correction at the
$0.05$ level. We see that the spatially-predicted test statistics increase
power with more detected changepoints. On the top right, we give the same plot
for Benjamini-Hochberg-corrected p-values based on spatial prediction. P-value
correction reduces the number of times changepoints were detected, yet more
changepoints are detected compared to only using data at each location.
Furthermore, the areas with changepoints detected for all five ensemble
correspond to the approximate latitude of the injection site ($15^{\circ}$N).
On the bottom right of Figure 7, we plot an estimate of
$\int_{0}^{1}\delta(\boldsymbol{s},u)du$, the averaged amount of change
associated with the epidemic time period. Grid cells were plotted if a
changepoint was detected using a Benjamini-Hochberg-corrected p-value. There
is generally an increase in stratospheric temperature across the ensembles
between 1-4 degrees K for detected locations in the tropics, in line with
expectations.

In general, the HSW++ data reinforces the validity of our approach to
spatially predict test statistics by detecting expected changes, and the form
of the changepoint model (AMOC or epidemic) can largely affect results and
interpretation of observed changes.

###  5.2 Application to MERRA-2 reanalysis

Figure 8: (Left) Longitudinal average of MERRA-2 reanalysis of stratospheric
temperature. (Top Right) Example of MERRA-2 reanalysis stratospheric
temperature data for $37.5^{\circ}\textrm{E}$ and $18^{\circ}\textrm{N}$.
(Bottom Right) Global average of stratospheric temperature.

We next turn to the MERRA-2 climate reanalysis data (GMAO,, 2015) which
assimilates historical observations with a climate model to provide data at
standardized locations and time points. We use stratospheric (near 50mb)
temperature data from 1984-1998 ($N=15$) to provide substantial data around
either side of the June 1991 eruption of Mt. Pinatubo. We use daily-averaged
data, wrangled into $m=365$ observations per year. In Figure 8, we plot an
overview of the MERRA-2 data. As expected, we see strong seasonal cycles that
dominates other sources of variability. In contrast to the HSW++ data, the
change in stratospheric temperature around the time of the aerosol injection
is less discernible. We use the same number of basis functions and principal
components as the HSW++ analysis.

Figure 9: MERRA-2 stratospheric temperature results for AMOC model. (Top Left)
Unadjusted p-values. (Top Right) Boxplot and points of the average change
estimate by $\hat{\tau}(\boldsymbol{s})$ for locations with significant BH-
adjusted p-value. (Bottom Left) Estimated $\hat{\tau}(\boldsymbol{s})$.
(Bottom Right) For each latitude, the proportion of locations (of 288
different longitudes) with $\hat{\tau}(\boldsymbol{s})$ for each year and
significant BH-corrected p-values.

We present results first from the AMOC model in Figure 9. There are low
p-values in the poles as well as around the equator. While detected
changepoints in the tropics (approximately $30^{\circ}$S to $30^{\circ}$N) may
be a direct result of the eruption of Mt. Pinatubo, changes near the poles
would be more likely associated with other changes or variability in the
climate. In the top right of Figure 9, we plot a boxplot of
$\hat{\Delta}(\boldsymbol{s}_{i})=\int_{0}^{1}\hat{\delta}(\boldsymbol{s}_{i},u)du$
for different $\boldsymbol{s}_{i}$ with a detected changepoint. The results
are mostly consistent with expected changes in stratospheric temperature based
on the eruption of Mt. Pinatubo. For example, for detected changes with the
change period beginning for 1992 or later ($\hat{\tau}(\boldsymbol{s})\geq
1991$), there is an expected decrease in stratospheric temperature on the
order of 1 to 2 K. This suggests the detection procedure detects a “return-to-
normal” type change in the data. Similarly, there are detected changepoints
with a later increase in stratospheric temperature in 1987 and 1990. Many
changes may also be associated with other sources of variability in the
climate data. In the bottom of Figure 9, we provide summaries of
$\hat{\tau}(\boldsymbol{s})$. Many changepoints are concentrated around the
time of the eruption, potentially detecting either a “return-to-normal” type
change or a change based on the original eruption.

Figure 10: (Left) Unadjusted epidemic p-values. (Right) Mean average change
estimate for the epidemic changepoint model for all locations by year of
detected change.

In Figure 10, we give results for the epidemic changepoint model. Since the
epidemic changepoint model is more complex, we have less power, and after
Benjamini-Hochberg correction to the p-values, no detections are made at the
$\alpha=0.05$ level. On the other hand, there are locations with significant
changes when not using p-value adjustment, and there are substantially lower
p-values in the tropics compared to the p-values from the AMOC change model.
As with the HSW++ data, the epidemic change model may be more appropriate. In
the right panel of Figure 10, we plot
$\textrm{mean}(\hat{\Delta}(\boldsymbol{s}))$ when grouping by both the
estimated $\hat{\tau}_{1}(\boldsymbol{s})$ and
$\hat{\tau}_{2}(\boldsymbol{s})$. The estimated changes again are in
accordance with the expected effects of Mt. Pinatubo on stratospheric
temperature. When $\hat{\tau}_{1}(\boldsymbol{s})\leq 1989$, an average
increase in stratospheric temperature on the order of $0.25-1.6$ K is
estimated, with these epidemic periods spanning the immediate aftermath of the
eruption (1991-1994). Alternately, when $\hat{\tau}_{1}(\boldsymbol{s})\geq
1991$, an average decrease in stratospheric temperature is estimated,
suggesting a return-to-normal-type changepoint is estimated using
$\hat{\tau}_{1}(\boldsymbol{s})$. More results for the epidemic model are
given in the supplement Section S5.

In Figure 11, we plot two examples of locations in our changepoint analysis.
We plot the MERRA-2 stratospheric temperature data and our estimates. Even
though the p-values are relatively large, the estimated changes use structure
in the data. On the left panel, the years 1989-1992, and especially the year
after the eruption of Mt. Pinatubo in June 1991, have a higher stratospheric
temperatures than other years. The AMOC changepoint model groups this 4-year
period with 1984-1988, leading to an estimated decrease in temperature after
the change in 1992. The epidemic changepoint selects 1989-1992 corresponding
to an expected increase in temperature. At the other location on the right
panel, both models group the higher temperatures immediately after the
eruption with 1984 through 1990, which had higher temperature compared to the
years 1993-1998. As with the HSW++ data, the changepoint detection procedures
may pick up unexpected events, and visualization and evaluation of the
estimated changes can inform resulting interpretation.

Figure 11: Two examples of stratospheric temperature and estimated change
points. The solid black line is the change detected using the AMOC model, with
p-values of $0.472$ (Left) and $0.996$ (Right). The two dashed red lines are
the change period detected using the epidemic model, with p-values of $0.428$
(Left) and $0.259$ (Right).

In the supplement Section S6, we also give the results on another atmospheric
variable from MERRA-2, aerosol optical depth (AOD). Overall, the results are
mostly similar to the results for stratospheric temperature. While power
continues to be a challenge for detecting changepoints, the estimated amount
of changes largely give expected changes associated with Mt. Pinatubo, while
describing complex spatial heterogeneity in the climate.

##  6 Discussion

In future efforts to understand changes in climate due to volcanic eruptions
or climate geoengineering events, the ability of statistical procedures to
detect and estimate such changes will play an important role in the
interpretation and response to such events. Such statistical procedures could
also play a role in the planning of climate interventions. We introduce an
approach using functional time series for such situations, and we
comprehensively evaluate its abilities and limitations for this purpose. In
this work, we provide a comprehensive framework for detecting changepoints in
globally-indexed functional time series that handles spatial heterogeneity in
the measured variable and its changepoints. We demonstrate that one can use
basic approaches to improve both changepoint detection power and the
estimation of the changepoint year using spatial information. We iteratively
apply the methodology on increasingly complex source data, demonstrating its
use on a statistical simulation study, a simplified climate model, and climate
reanalysis data. Moreover, we use the simulation study to evaluate different
test statistics and approaches for using spatial information, finding that
estimating the model under the null hypothesis and using score-based test
statistics consistently improve estimates.

The full-variability analysis on MERRA-2 reanalysis data is challenging due to
two main reasons. First, there are complex sources of variability. While the
functional time series approach handles seasonality naturally, other
variability like long-term trends are ignored. The second challenge is a lack
of power due to a reduced sample size. In our analysis on climate reanalysis
data, we use $N=15$ years; in this setting, it is challenging to detect
changes based on asymptotic results that take $N\to\infty$. However, this
problem does not have a simple solution. One could aim to increase the number
of years of data, yet over longer periods of time different events may be
influential. For example, the 1982 volcanic eruption of El Chinchón in Mexico
likely influenced the climate in similar ways to the eruption of Mt. Pinatubo
in 1991 (Marshall et al.,, 2022). One could also aim to use time periods more
granular than years for the functional time series, yet one would need an
alternate approach to handle seasonality. One could also only infer
changepoints at a reduced set of locations, reducing the challenge of the
multiple testing problem, or apply more sophisticated spatial testing
approaches (Risser et al.,, 2019; Yun et al.,, 2022). This challenge is
especially pronounced for the epidemic changepoint model, as the epidemic
change period may only last one or two years. Therefore, considering the
problem under the framework of anomaly detection may be useful. To evaluate
the test statistics’ sensitivity to the reduced sample size, we repeated the
HSW++ analysis for the first ensemble member using $N=12$ and $m=400$ instead
of $N=40$ and $m=120$, with results presented in the supplement Section S4. As
expected, this decreased the power for detecting changepoints, suggesting that
the reduced $N$ substantially challenged the full-variability analysis.

One other limitation of this work is the lack of dependence assumed across
$k$. As detailed in Lund et al., (2022), ignoring such dependence may lead to
conservative or anticonservative testing results. However, since, the
functional time series are partitioned into relatively long time periods of
years, dependence may be weaker between years compared to other settings
(where $k$ might index weeks or days). As written, the spatial approaches
developed here should be especially useful when the variable varies relatively
smoothly without hard boundaries between regions of change and no change.
Otherwise, the spatial model may unduly extrapolate changes in the mean to
areas where there is no true change in the mean. Using spatial approaches like
clustering would be useful for this setting with a variable like surface
temperature, which is strongly influenced by hard topographic boundaries.

## Acknowledgment

This paper describes objective technical results and analysis. Any subjective
views or opinions that might be expressed in the paper do not necessarily
represent the views of the U.S. Department of Energy or the United States
Government. This work was supported by the Laboratory Directed Research and
Development program at Sandia National Laboratories; a multi-mission
laboratory managed and operated by National Technology and Engineering
Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell
International, Inc., for the U.S. Department of Energy’s National Nuclear
Security Administration under contract DE-NA0003525.

## References

  * Alegría et al., (2020) Alegría, A., Emery, X., & Lantuéjoul, C. (2020).  The turning arcs: a computationally efficient algorithm to simulate isotropic vector-valued gaussian random fields on the d-sphere.  Statistics and Computing, 30, 1403–1418.
  * Aston & Kirch, (2012) Aston, J. A. & Kirch, C. (2012).  Detecting and estimating changes in dependent functional data.  Journal of Multivariate Analysis, 109, 204–220.
  * Aue et al., (2009) Aue, A., Gabrys, R., Horváth, L., & Kokoszka, P. (2009).  Estimation of a change-point in the mean function of functional data.  Journal of Multivariate Analysis, 100(10), 2254–2269.
  * Aue et al., (2018) Aue, A., Rice, G., & Sönmez, O. (2018).  Detecting and dating structural breaks in functional data without dimension reduction.  JRSSB, 80(3), 509–529.
  * Aue et al., (2020) Aue, A., Rice, G., & Sönmez, O. (2020).  Structural break analysis for spectrum and trace of covariance operators.  Environmetrics, 31(1), e2617.  e2617 env.2617.
  * Benjamini & Heller, (2007) Benjamini, Y. & Heller, R. (2007).  False discovery rates for spatial signals.  JASA, 102(480), 1272–1281.
  * Benjamini & Hochberg, (1995) Benjamini, Y. & Hochberg, Y. (1995).  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  JRSSB, 57(1), 289–300.
  * Berkes et al., (2009) Berkes, I., Gabrys, R., Horvath, L., & Kokoszka, P. (2009).  Detecting changes in the mean of functional observations.  JRSSB, 71, 927–946.
  * Bluth et al., (1992) Bluth, G. J. S., Doiron, S. D., Schnetzler, C. C., Krueger, A. J., & Walter, L. S. (1992).  Global tracking of the SO2 clouds from the June, 1991 Mount Pinatubo eruptions.  Geophysical Research Letters, 19(2), 151–154.
  * Cao et al., (2022) Cao, J., Zhang, J., Sun, Z., & Katzfuss, M. (2022).  Locally anisotropic covariance functions on the sphere.
  * Caponera et al., (2022) Caponera, A., Fageot, J., Simeoni, M., & Panaretos, V. M. (2022).  Functional estimation of anisotropic covariance and autocovariance operators on the sphere.  Electronic Journal of Statistics, 16(2), 5080 – 5148.
  * Chiou et al., (2019) Chiou, J.-M., Chen, Y.-T., & Hsing, T. (2019).  Identifying multiple changes for a functional data sequence with application to freeway traffic segmentation.  The Annals of Applied Statistics, 13(3), 1430 – 1463.
  * Dai & Müller, (2018) Dai, X. & Müller, H.-G. (2018).  Principal component analysis for functional data on Riemannian manifolds and spheres.  The Annals of Statistics, 46(6B), 3334 – 3361.
  * Dette & Kutta, (2021) Dette, H. & Kutta, T. (2021).  Detecting structural breaks in eigensystems of functional time series.  Electronic Journal of Statistics, 15(1), 944 – 983.
  * Dette & Quanz, (2022) Dette, H. & Quanz, P. (2022).  Detecting relevant changes in the spatiotemporal mean function.  Journal of Time Series Analysis, n/a(n/a).
  * García-Portugués & Verdebout, (2018) García-Portugués, E. & Verdebout, T. (2018).  An overview of uniformity tests on the hypersphere.  arXiv:1804.00286.
  * García-Portugués & Verdebout, (2021) García-Portugués, E. & Verdebout, T. (2021).  sphunif: Uniformity Tests on the Circle, Sphere, and Hypersphere.  R package version 1.0.1.
  * GMAO, (2015) GMAO (2015).  MERRA-2 tavg3_3d_asm_nv: 3d, 3-hourly, time-averaged, model-level, assimilation, assimilated meteorological fields v5.12.4.
  * Gromenko et al., (2017) Gromenko, O., Kokoszka, P., & Reimherr, M. (2017).  Detection of change in the spatiotemporal mean function.  JRSSB, 79(1), 29–50.
  * Harris et al., (2022) Harris, T., Li, B., & Tucker, J. D. (2022).  Scalable multiple changepoint detection for functional data sequences.  Environmetrics, 33(2), e2710.
  * Held & Suarez, (1994) Held, I. & Suarez, M. (1994).  A proposal for the intercomparison of the dynamical cores of atmospheric general circulation models.  Bulletin of the American Meteorological Society, 75(10), 1825–1830.
  * Hollowed, (2022) Hollowed, J. (2022).  A simplified sub-grid parameterization for volcanic aerosol injections of the stratosphere in E3SMv2 in support of the DOE CLDERA project.  Preliminary exam research prospectus, (pp. 1–17).  https://jhollowed.github.io/pages/content/assets/docs/hollowed_prelim.pdf.
  * Hörmann & Kokoszka, (2012) Hörmann, S. & Kokoszka, P. (2012).  Functional time series.  In Handbook of statistics, volume 30 (pp. 157–186). Elsevier.
  * Horváth et al., (2022) Horváth, L., Rice, G., & Zhao, Y. (2022).  Change point analysis of covariance functions: A weighted cumulative sum approach.  Journal of Multivariate Analysis, 189, 104877.
  * Jeong et al., (2017) Jeong, J., Jun, M., & Genton, M. G. (2017).  Spherical process models for global spatial statistics.  Statistical Science, 32(4), 501–513, 13.
  * (26) Jiao, S., Chan, N.-H., & Yau, C.-Y. (2022a).  Enhanced change-point detection in functional means.
  * (27) Jiao, S., Frostig, R. D., & Ombao, H. (2022b).  Break point detection for functional covariance.  Scandinavian Journal of Statistics, n/a(n/a).
  * Katzfuss et al., (2020) Katzfuss, M., Jurek, M., Zilber, D., & Gong, W. (2020).  GPvecchia: Scalable Gaussian-Process Approximations.  R package version 0.1.3.
  * Khodri et al., (2017) Khodri, M., Izumo, T., Vialard, J., Janicot, S., Cassou, C., Lengaigne, M., Mignot, J., Gastineau, G., Guilyardi, E., Lebas, N., et al. (2017).  Tropical explosive volcanic eruptions can trigger El Niño by cooling tropical Africa.  Nature communications, 8(1), 778.
  * Krivoruchko & Gribov, (2020) Krivoruchko, K. & Gribov, A. (2020).  Distance metrics for data interpolation over large areas on earth’s surface.  Spatial Statistics, 35, 100396.
  * Labitzke & McCormick, (1992) Labitzke, K. & McCormick, M. (1992).  Stratospheric temperature increases due to Pinatubo aerosols.  Geophysical Research Letters, 19(2), 207–210.
  * Lenton & Vaughan, (2009) Lenton, T. M. & Vaughan, N. E. (2009).  The radiative forcing potential of different climate geoengineering options.  Atmospheric Chemistry and Physics, 9(15), 5539–5561.
  * Lund et al., (2022) Lund, R. B., Beaulieu, C., Killick, R., Lu, Q., & Shi, X. (2022).  Good practices and common pitfalls in climate time series changepoint techniques: a review.
  * Majumdar et al., (2005) Majumdar, A., Gelfand, A. E., & Banerjee, S. (2005).  Spatio-temporal change-point modeling.  Journal of Statistical Planning and Inference, 130(1-2), 149–166.
  * Marshall et al., (2022) Marshall, L. R., Maters, E. C., Schmidt, A., Timmreck, C., Robock, A., & Toohey, M. (2022).  Volcanic effects on climate: recent advances and future avenues.  Bulletin of Volcanology, 84(5), 54.
  * Moradi et al., (2023) Moradi, M., Cronie, O., Pérez-Goya, U., & Mateu, J. (2023).  Hierarchical spatio-temporal change-point detection.  The American Statistician, 0(0), 1–11.
  * NASEM, (2021) NASEM (2021).  Reflecting sunlight: recommendations for solar geoengineering research and research governance.  The National Academies Press.
  * Paciorek & Schervish, (2006) Paciorek, C. J. & Schervish, M. J. (2006).  Spatial modelling using a new class of nonstationary covariance functions.  Environmetrics: The official journal of the International Environmetrics Society, 17(5), 483–506.
  * Paik & Min, (2018) Paik, S. & Min, S.-K. (2018).  Assessing the impact of volcanic eruptions on climate extremes using CMIP5 models.  Journal of Climate, 31(14), 5333 – 5349.
  * Parker et al., (1996) Parker, D., Wilson, H., Jones, P. D., Christy, J., & Folland, C. K. (1996).  The impact of mount pinatubo on world-wide temperatures.  International Journal of Climatology: A Journal of the Royal Meteorological Society, 16(5), 487–497.
  * Porcu et al., (2018) Porcu, E., Alegria, A., & Furrer, R. (2018).  Modeling temporally evolving and spatially globally dependent data.  International Statistical Review, 86(2), 344–377.
  * Porcu et al., (2019) Porcu, E., Castruccio, S., Alegria, A., & Crippa, P. (2019).  Axially symmetric models for global data: a journey between geostatistics and stochastic generators.  Environmetrics, 30(1), e2555.
  * Reeves et al., (2007) Reeves, J., Chen, J., Wang, X. L., Lund, R., & Lu, Q. Q. (2007).  A review and comparison of changepoint detection techniques for climate data.  Journal of Applied Meteorology and Climatology, 46(6), 900 – 915\.
  * Rice & Zhang, (2022) Rice, G. & Zhang, C. (2022).  Consistency of binary segmentation for multiple change-point estimation with functional data.  Statistics & Probability Letters, 180, 109228.
  * Risser et al., (2019) Risser, M. D., Paciorek, C. J., & Stone, D. A. (2019).  Spatially dependent multiple testing under model misspecification, with application to detection of anthropogenic influence on extreme climate events.  JASA, 114(525), 61–78.
  * Stein, (2005) Stein, M. L. (2005).  Nonstationary spatial covariance functions.  Unpublished technical report.
  * Stoehr et al., (2021) Stoehr, C., Aston, J. A. D., & Kirch, C. (2021).  Detecting changes in the covariance structure of functional time series with application to fMRI data.  Econometrics and Statistics, 18, 44–62.
  * Su et al., (2014) Su, J., Kurtek, S., Klassen, E., & Srivastava, A. (2014).  Statistical analysis of trajectories on Riemannian manifolds: bird migration, hurricane tracking and video surveillance.
  * Tucker & Yarger, (2022) Tucker, J. D. & Yarger, D. (2022).  Elastic functional changepoint detection.  https://arxiv.org/abs/2211.12687.
  * Wahba, (1990) Wahba, G. (1990).  Spline models for observational data.  SIAM.
  * Wang et al., (2022) Wang, M., Harris, T., & Li, B. (2022).  Asynchronous changepoint estimation for spatially correlated functional time series.  Journal of Agricultural, Biological and Environmental Statistics.
  * Wegner & Wendler, (2022) Wegner, L. & Wendler, M. (2022).  Robust change-point detection for functional time series based on U-statistics and dependent wild bootstrap.
  * Williamson et al., (1998) Williamson, D. L., Olson, J. G., & Boville, B. A. (1998).  A comparison of semi-Lagrangian and Eulerian tropical climate simulations.  Monthly Weather Review, 126(4), 1001 – 1012.
  * Yun et al., (2022) Yun, S., Zhang, X., & Li, B. (2022).  Detection of local differences in spatial characteristics between two spatiotemporal random fields.  JASA, 117(537), 291–306.
  * Zhao et al., (2019) Zhao, Z., Ma, T. F., Ng, W. L., & Yau, C. Y. (2019).  A composite likelihood-based approach for change-point detection in spatio-temporal process.  arXiv preprint arXiv:1904.06340.