manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
Parameter uncertainty quantification in an idealized GCM with
1
a seasonal cycle
2
Michael F. Howland
1
,
2
, Oliver R. A. Dunbar
2
, Tapio Schneider
2
3
1
Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, MA, USA
4
2
California Institute of Technology, Pasadena, CA, USA
5
Corresponding author: Michael F. Howland,
mhowland@mit.edu
–1–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
Abstract
6
Climate models are generally calibrated manually by comparing selected climate statistics, such
7
as the global top-of-atmosphere energy balance, to observations. The manual tuning only tar-
8
gets a limited subset of observational data and parameters. Bayesian calibration can estimate
9
climate model parameters and their uncertainty using a larger fraction of the available data and
10
automatically exploring the parameter space more broadly. In Bayesian learning, it is natu-
11
ral to exploit the seasonal cycle, which has large amplitude, compared with anthropogenic cli-
12
mate change, in many climate statistics. In this study, we develop methods for the calibration
13
and uncertainty quantification (UQ) of model parameters exploiting the seasonal cycle, and
14
we demonstrate a proof-of-concept with an idealized general circulation model (GCM). Un-
15
certainty quantification is performed using the calibrate-emulate-sample approach, which com-
16
bines stochastic optimization and machine learning emulation to speed up Bayesian learning.
17
The methods are demonstrated in a perfect-model setting through the calibration and UQ of
18
a convective parameterization in an idealized GCM with a seasonal cycle. Calibration and UQ
19
based on seasonally averaged climate statistics, compared to annually averaged, reduces the
20
calibration error by up to an order of magnitude and narrows the spread of posterior distri-
21
butions by factors between two and five, depending on the variables used for UQ. The reduc-
22
tion in the size of the parameter posterior distributions leads to a reduction in the uncertainty
23
of climate model predictions.
24
Plain Language Summary
25
Climate models rely on parameterizations of physical processes that cannot be resolved
26
with available computational resources. Empirical representations of physical processes, such
27
as turbulence and cloud physics, reduce the computational cost of simulations, but introduce
28
new unknown parameters into the climate model. The unknown parameters contribute to un-
29
certainties associated with climate model predictions. Historically, fixed values of the model
30
parameters have been hand-tuned using scientific intuition and a limited amount of available
31
data. We develop methods for the computationally efficient estimation of the unknown climate
32
model parameters and their uncertainty systematically from data, by using gradient-free op-
33
timization and machine learning. Many processes and observable statistics of Earth’s climate
34
used to produce this data are influenced by the seasonal cycle. We demonstrate that the in-
35
corporation of seasonal information into these statistics significantly improves the resulting cal-
36
ibration of climate model parameters, in contrast to annually averaged information alone. We
37
show that including seasonal information also reduces the uncertainty associated with the model
38
parameters, which consequently reduces the uncertainty of climate model predictions.
39
1 Introduction
40
The objective of quantifying uncertainty in computational models arises in a wide range
41
of applications, including weather and climate modeling (Schneider, Lan, et al., 2017), fluid
42
dynamics (Duraisamy et al., 2019), and energy systems (Constantinescu et al., 2010). Often,
43
uncertainty associated with predictions from computational models is the result of processes
44
that cannot be resolved on the computational grid, either due to computational complexity lim-
45
itations (Meneveau & Katz, 2000) or due to uncertainty associated with the process itself (Schneider,
46
Teixeira, et al., 2017). In general circulation models (GCMs), primary uncertainties arise from
47
the representation of subgrid-scale turbulence, convection, and cloud physics, which have a
48
significant impact on the evolution of climate under rising greenhouse gases (Cess et al., 1989;
49
Bony & Dufresne, 2005; Webb et al., 2013; Suzuki et al., 2013; Brient & Schneider, 2016).
50
While clouds are associated with turbulence and convective updrafts with scales of
O
(10 m)
,
51
modern climate simulations have a typical horizontal resolution of
O
(10 km)
–
O
(100 km)
(Delworth
52
et al., 2012; Hoegh-Guldberg et al., 2018). Climate simulations rely on physically motivated
53
parameterizations that model the effects of subgrid-scale processes such as clouds and turbu-
54
lence on the resolved scales (Hourdin et al., 2013). Such parameterizations come with para-
55
–2–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
metric and structural uncertainties; quantifying these uncertainties and how they percolate into
56
climate projections remains an outstanding challenge (Schneider, Teixeira, et al., 2017).
57
Physical parameterizations have historically been individually developed and calibrated
58
using data from isolated experiments (Smagorinsky, 1963; Spalart & Allmaras, 1992; Jakob,
59
2010; Golaz et al., 2013; Hourdin et al., 2017). They are further adjusted so that global mod-
60
els that incorporate them satisfy selected large-scale observational or physical constraints, such
61
as a closed top-of-atmosphere energy balance or reproduction of the 20th-century global-mean
62
temperature evolution (Mauritsen et al., 2012; Hourdin et al., 2017; Schmidt et al., 2017). Model
63
calibration is usually done manually, focusing on a small subset of model parameters and ex-
64
ploiting only a fraction of the available observational data.
65
As a step toward automating and augmenting this process, here we further develop al-
66
gorithms for model calibration and uncertainty quantification (UQ) that in principle allow mod-
67
els to learn from large datasets and that scale to high-dimensional parameter spaces. In pre-
68
vious work, these algorithms have been demonstrated in simple conceptual models (Cleary et
69
al., 2021) and in a statistically stationary idealized GCM (Dunbar et al., 2020). We take the
70
next step and demonstrate how these algorithms can exploit seasonal variations, which for many
71
climate statistics are large relative to the climate changes expected in the coming decades and
72
contain exploitable information about the response of the climate system to perturbations (Knutti
73
et al., 2006; Schneider et al., 2021).
74
Whereas numerical weather prediction assimilates spatiotemporally evolving trajectories
75
of atmospheric states as initial conditions for forecasts (e.g., Kalnay (2003); Bauer et al. (2015)),
76
in climate modeling it is preferable to assimilate time-averaged climate statistics. This focuses
77
the learning problem on quantities of interest in climate predictions (i.e., climate statistics, in-
78
cluding higher-order statistics such as precipitation extremes), and it avoids the need to esti-
79
mate uncertain atmospheric initial conditions on which trajectories of states depend (Cleary
80
et al., 2021). Calibration and UQ of climate models on the basis of time-averaged statistics
81
smooths the prediction-error based objective function and enables the use of data that have dif-
82
ferent resolution than the climate simulations (Schneider, Lan, et al., 2017; Dunbar et al., 2020).
83
Given the dual desires to avoid having to estimate atmospheric initial conditions, which are
84
forgotten over about 2 weeks (Zhang et al., 2019), and to exploit seasonal variations, it be-
85
comes natural to choose averaging timescales between around 30 and 90 days. Such averag-
86
ing timescales are the focus of this study.
87
We consider the calibration and UQ of convective parameters in an idealized GCM with
88
seasonally varying insolation (Frierson et al., 2006; O’Gorman & Schneider, 2008; Bordoni
89
& Schneider, 2008). We develop an extension of the calibrate-emulate-sample (CES) Bayesian
90
learning methodology (Cleary et al., 2021) to enable the use of statistics computed from a non-
91
stationary statistical state. The qualitative and quantitative impacts of the time-averaging length
92
of the climate statistics on the parameter calibration and UQ are assessed. GCMs are gener-
93
ally tuned in situations which have low parameter identifiability based on available climate data
94
(Hourdin et al., 2017). We perform numerical experiments with observable climate statistics
95
that are informative about the convective parameters we wish to calibrate, but are not in any
96
simple and direct way related to them. We also explore less informative climate statistics, which
97
highlight the benefits of incorporating seasonal variations in the climate statistics that are be-
98
ing exploited for UQ.
99
The remainder of this paper is organized as follows. In section 2, the Bayesian learn-
100
ing methods for time-dependent problems are introduced. The numerical details of the sea-
101
sonally forced GCM and UQ experiments are introduced in section 3. The calibration and UQ
102
results are shown in section 4, and conclusions are provided in section 5.
103
–3–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
2 Uncertainty quantification methods
104
The goal of this study is to estimate the probability distributions associated with model
105
parameters
θ
that are used by a GCM and about which only imprecise prior information is
106
known. We will consider a Bayesian approach to the estimation of the probability distributions
107
of model parameters
θ
, where we seek
P
(
θ
|
y
)
, the conditional probability distribution of pa-
108
rameters
θ
given observed data
y
. The GCM is a computationally expensive numerical model
109
that evolves climate states in time. From the GCM, we extract statistical information that is
110
denoted by
G
(
θ
)
. Here,
G
(
θ
)
includes the numerical integration of the GCM states and the
111
aggregation of relevant climate statistics in time. Although the data
y
will, in general, not pro-
112
vide direct information about
θ
, we can estimate
P
(
θ
|
y
)
using Bayesian learning by compar-
113
ing the GCM outputs
G
(
θ
)
with data
y
.
114
Standard methods for Bayesian learning include Markov Chain Monte Carlo (MCMC)
115
(Brooks et al., 2011), which typically requires
O
(10
5
)
forward model evaluations of
G
(
θ
)
to
116
sample the posterior distribution (Geyer, 2011). Instead, we perform calibration and UQ us-
117
ing the recently developed CES methodology (Cleary et al., 2021), which consists of three steps:
118
(1) Ensemble Kalman processes (Schillings & Stuart, 2017; Garbuno-Inigo et al., 2020) are
119
used to calibrate parameters and to generate input-output pairs of the mapping
θ
7→ G
(
θ
)
;
120
(2) Gaussian process (GP) regression is used to train an emulator
G
GP
(
θ
)
of the mapping
θ
7→
121
G
(
θ
)
using the training points generated in the calibration step; and (3) MCMC sampling with
122
the computationally efficient GP emulator
G
GP
(
θ
)
rather than the expensive forward model
123
G
(
θ
)
is used to estimate the posterior distribution
P
(
θ
|
y
)
. The CES methodology has previ-
124
ously been used for calibration and UQ of parameters in simple model problems such as Darcy
125
flow and Lorenz systems (Cleary et al., 2021) and for convective parameters in a statistically
126
stationary GCM (Dunbar et al., 2020). More recent methodological developments by Lan et
127
al. (2021) enabled the CES framework to perform simultaneous UQ on
O
(1000)
parameters
128
using deep neural network-based emulation and MCMC suited to high-dimensional spaces (Beskos
129
et al., 2008, 2011). A schematic of the CES methodology is shown in Figure 1.
130
The Bayesian learning methodology used in this study is described in the following sec-
131
tions. In section 2.1, the inverse problem of estimating
P
(
θ
|
y
)
in a setting with a periodic cy-
132
cle is introduced. In section 2.2, the ensemble Kalman process calibration method is outlined.
133
Section 2.3 introduces the GP emulation in an uncorrelated transformed space, obtained by
134
a singular value decomposition (principal component analysis) of the noise covariance matrix.
135
Section 2.4 describes how the Bayesian posterior distribution is approximated using the GP
136
emulator.
137
2.1 Seasonal GCM inverse problem
138
With fixed insolation (Frierson et al., 2006; O’Gorman & Schneider, 2008), GCM statis-
tics are stationary and ergodic. With the seasonal cycle incorporated, the insolation varies as
a function of the ordinal day in the simulation. The resulting GCM states are statistically cy-
clostationary, with a dependence on the ordinal day. The data we use are constructed account-
ing for the seasonally varying boundary conditions, as the time-averaged GCM statistics
G
T
(
θ
,ξ,t
) =
1
T
∫
t
+
T
t
H
(
t
′
;
θ
,ξ
)
dt
′
.
(1)
Here,
H
(
t
;
θ
,ξ
)
∈
R
N
y
represents the output states of the GCM, depending on time
t
. The
139
number of states being measured is
N
y
. The integration length is specified as
T
. The initial
140
conditions for the GCM integration are represented by
ξ
, which depend on the initial time
t
.
141
The time-averaged statistics are
G
T
(
θ
,ξ,t
)
∈
R
N
y
.
142
With seasonally varying insolation, the time-averaged statistics in Eq. (1) depend on the
ordinal day (
t
) in the GCM. With a specified integration length
T
, there are a corresponding
–4–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
Calibrate
EKI/EKS using
Transformation
{
θ
i
,
G
(
θ
i
)
}
i
→
{
θ
i
,
̃
G
(
θ
i
)
}
i
Emulate
Gaussian process
{
̃
G
(
θ
i
)
,
̃
Σ
(
θ
i
)
}
i
→
{
̃
G
GP
(
θ
i
)
,
̃
Σ
GP
(
θ
i
)
}
i
{
θ
i
,
G
(
θ
i
)
}
i
Uncorrelated space from truncated SVD
Sample
MCMC using
Φ
MCMC
(
θ
,
̃
y
)
̃
G
GP
P
(
θ
|
̃
y
)
Calibrate
−
Emulate
−
Sample
Φ
(
θ
,
y
)
Inputs
:
•
Data
y
, covariance
Σ
•
Characteristic values
y
c
•
Prior
P
(
θ
)
Output
:
̃
Σ
GP
Figure 1: Schematic of the calibrate-emulate-sample (CES) methodology to estimate model pa-
rameters
θ
. With inputs of data
y
, noise covariance
Σ
, characteristic values of the data
y
c
(for
normalization), and the prior
P
(
θ
)
, the calibration stage generates input-output pairs
{
θ
i
,
G
(
θ
i
)
}
i
.
The input-output mapping is emulated using Gaussian process regression in a transformed, uncor-
related space (
̃
·
), obtained from a truncated singular value decomposition on the noise covariance
matrix
Σ
. The GP emulator is used for efficient sampling using MCMC to approximate the pos-
terior distribution
P
(
θ
|
̃
y
)
. The objective functions for calibration and sampling are denoted by
Φ(
θ
,
y
)
and
Φ
MCMC
(
θ
,
̃
y
)
, respectively.
set of time-averaged statistics
G
T,j
(
θ
,ξ
) =
1
T
∫
t
0
+
jT
t
0
+(
j
−
1)
T
H
(
t
′
;
θ
,ξ
)
dt
′
,
(2)
where
j
= 1
through
N
s
is an index representing the time of year in the GCM simulation,
143
starting from
j
= 1
at vernal equinox (time
t
0
). The integration windows are non-overlapping.
144
The length of the year in the GCM is
T
yr
= 360 d
, and the resulting number of non-overlapping
145
time-averaged statistics are
N
s
=
T
yr
/T
. In this framework, the number of batches of statis-
146
tics
N
s
is a design parameter set by the integration timescale
T
. For
T
= 90 d
, the data are
147
aggregated seasonally, with
N
s
= 4
. For
T
= 360 d
, the statistics are averaged over the
148
full year in the GCM, corresponding to annually averaged climate statistics, and
N
s
= 1
.
149
The GCM data are the concatenation of the time-averaged batches over collated ordi-
nal days
G
(
θ
,ξ
) = [
G
T,
1
(
θ
,ξ
1
)
,
G
T,
2
(
θ
,ξ
2
)
,...,
G
T,N
s
(
θ
,ξ
N
s
)]
.
(3)
In this study, we will consider perfect-model numerical experiments, with the data
y
being con-
structed using the same GCM at target true parameters
θ
=
θ
†
. The synthetic data are given
by
y
= [
G
T,
1
(
θ
†
,ξ
1
)
,
G
T,
2
(
θ
†
,ξ
2
)
,...,
G
T,N
s
(
θ
†
,ξ
N
s
)]
.
(4)
The size of
y
and
G
is
R
N
, where
N
=
N
s
·
N
y
. The resulting inverse problem relating pa-
rameters and data is
y
=
G
(
θ
,ξ
) +
η,
(5)
where
η
∼
N
(0
,
∆)
is a realization of normal measurement error with zero mean and co-
150
variance matrix
∆
. We generate synthetic data and forward model data starting from the same
151
initial conditions
ξ
in Eq. (5). Since the states of the GCM depend on the ordinal day, the av-
152
eraging operation in Eq. (2) must be consistently aligned in ordinal day between the synthetic
153
data and GCM outputs.
154
The states of the GCM depend on the boundary conditions, specific to the ordinal day,
and the initial conditions
ξ
. The ensemble average of independent realizations of
G
(
θ
,ξ
i
)
over
–5–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
different initial conditions
ξ
i
is
G
∞
(
θ
) = lim
M
→∞
1
M
M
∑
i
=1
G
(
θ
,ξ
i
)
.
(6)
In the limit of infinite realizations of the climate statistics, the dependence of the ensemble
155
average
G
∞
(
θ
)
on the initial conditions trends to zero by the central limit theorem. The cen-
156
tral limit theorem applies since
G
(
θ
,ξ
i
)
contains a full year of data by construction (Eq. 3)
157
and is therefore a statistically stationary object with respect to the cycle. But due to compu-
158
tational limitations, only a finite number of ensemble realizations is available in practice. For
159
the synthetic data, we can similarly define
y
∞
; however,
y
∞
is generally not accessible from
160
observations.
161
The atmospheric initial condition is of minimal practical value in the climate model set-
ting, and we aim to avoid its estimation. We therefore reformulate the inverse problem such
that
G
and
y
can have different initial conditions, and such that we do not require the estima-
tion of the initial conditions. To do so, we assume the forward model output
G
(
θ
,ξ
)
, is a noisy,
finite average approximation of the infinite ensemble average of the climate statistics, that is,
G
(
θ
,ξ
) =
G
∞
(
θ
) +
N
(0
,
Σ)
, where
Σ
is the internal variability covariance matrix of the
GCM (Cleary et al., 2021). The realizations
G
(
θ
,ξ
)
and
y
are both subject to the internal vari-
ability of the climate system. In this study, given the definition of
y
and
G
(
θ
,ξ
)
, the internal
variability is the interannual variability of the GCM. As a result, the inverse problem becomes
y
=
G
∞
(
θ
) +
γ,
(7)
where
γ
∼
N
(0
,
∆ + Σ)
. In Eq. (7), the initial conditions are removed from the inference
problem. The synthetic data are collected into the matrix
Y
∈
R
N
×
N
d
, where
N
is the di-
mension of the data space and
N
d
is the number of data samples, or years in this study. The
internal variability covariance matrix is computed from the synthetic data matrix
Y
such that
Σ = cov(
Y
)
. We assume that
Σ
does not vary as a function
θ
such that
Σ(
θ
†
)
≈
Σ(
θ
)
and
that
Σ
does not depend on initial conditions
ξ
. Invoking Gaussian error assumptions based on
the central limit theorem, the corresponding negative log-likelihood objective function then is
Φ(
θ
,
y
) =
1
2
‖
y
−G
∞
(
θ
)
‖
2
∆+Σ
(8)
where
‖·‖
A
=
‖
A
−
1
/
2
·‖
2
. The likelihood is (Kaipio & Somersalo, 2006)
P
(
θ
|
y
)
∝
exp(
−
Φ(
θ
,
y
))
.
(9)
In this perfect-model setting, where the synthetic data and forward model are obtained
from the same GCM, there is no direct measurement error (
η
) and no systematic model er-
ror. To emulate measurement error, we add Gaussian noise to the GCM output, with zero mean
and covariance matrix
∆ = diag(
δ
2
k
)
.
(10)
Here, the noise standard deviation is defined such that
δ
k
= min
(
C
min
[
dist(
y
k
+ 2
√
Σ
kk
,∂
Ω
k
)
,
dist(
y
k
−
2
√
Σ
kk
,∂
Ω
k
)
]
,C
m
·
y
k
)
,
(11)
where
C
m
= 0
.
1
caps the maximum measurement error standard deviation to
10%
of the mean
162
data values and
C
= 0
.
2
controls how close the noise-perturbed data can come within phys-
163
ical boundaries
∂
Ω
i
of the data (e.g., to keep relative humidities between 0 and 1 with high
164
probability).
165
2.2 Calibrate: Ensemble Kalman Inversion
166
The first stage of the CES methodology is to calibrate the parameters
θ
of the model
167
based on data
y
. We perform calibration with independent realizations of
G
(
θ
i
,ξ
i
)
, viewed
168
–6–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
as noisy approximations of
G
∞
(
θ
i
)
. Calibration is performed using ensemble Kalman meth-
169
ods, which demonstrate theoretical success, in idealized problems, and empirical success, in
170
complex problems, to optimize parameters under such noise (Duncan et al., 2021).
171
The utility of the calibration stage is two-fold: (1) optimize parameters to minimize the
mismatch between model output and data; and (2) provide good parameter–model output pairs
(
θ
i
,
G
(
θ
i
,ξ
i
))
for training an emulator of the parameter-to-data map, with a higher density of
training points near the optimal parameters. The ensemble Kalman filter (EnKF) is a Kalman
filter implementation in which the covariances are approximated using Monte Carlo sampling
(Evensen, 2003). The EnKF has been used widely for derivative-free state estimation in nu-
merical weather prediction (e.g., Houtekamer and Zhang (2016)) and model-based control (e.g.,
Howland, Ghate, et al. (2020)). Ensemble Kalman methods for Bayesian inversion were in-
troduced by Chen and Oliver (2012) and Emerick and Reynolds (2013). These methods prov-
ably draw samples from the posterior of linear inverse problems subject to additive Gaussian
noise; however, they fail to do so more in more general, nonlinear problems. Recognizing this,
the offline ensemble Kalman inversion (EKI) (Iglesias et al., 2013) algorithm was introduced
for classical, optimization-based inversion. EKI generally drives the ensemble members toward
consensus near the optimal solution of the inverse problem (Schillings & Stuart, 2017). The
parameter update of ensemble member
m
at iteration step
n
is
θ
(
n
+1)
m
=
θ
(
n
)
m
+
C
(
n
)
θ
G
(
Σ + ∆ +
C
(
n
)
GG
)
−
1
(
y
−G
(
θ
(
n
)
m
,ξ
m
)
)
,
(12)
where
C
θ
G
is the empirical cross-covariance between the parameters and the model outputs,
172
and
C
GG
is the empirical covariance of the model outputs.
173
EKI is guaranteed to find the optimal parameters in linear problems (Schillings & Stu-
174
art, 2017) and has empirical success in nonlinear problems (Dunbar et al., 2020). Several other
175
approaches exist for estimation in nonlinear problems. Li and Reynolds (2009) and Sakov et
176
al. (2012) developed iterative EnKF methods which improve state estimation in strongly non-
177
linear problems. However, the distribution of the ensemble does not converge to the posterior
178
distribution in the limit of infinite members for nonlinear problems (Zhou et al., 2006; An-
179
nan & Hargreaves, 2007; Le Gland et al., 2009; Garbuno-Inigo et al., 2020), necessitating the
180
emulation and sampling described in the following sections.
181
The number of ensemble members and EKI iterations are hyperparameters; they are set
182
to standard values of
N
ens
= 100
and
N
it
= 5
in this study (Schillings & Stuart, 2017; Dun-
183
bar et al., 2020; Cleary et al., 2021). Each ensemble member is run through the GCM, ini-
184
tialized from the same initial conditions. A spin-up period of one year (
360
days) is run be-
185
fore the statistics are computed using Eq. (2), to ensure the forward model realizations are sub-
186
ject to differing instantiations of internal variability in the chaotic GCM system.
187
2.3 Emulate: Gaussian process emulation
188
We emulate the mapping from parameters to model output using a machine learning method
189
that enables the rapid execution of the mapping, compared to the computationally expensive
190
forward model. The calibration stage (section 2.2) results in
N
t
=
N
ens
·
N
it
input-output
191
pairs
{
θ
i
,
G
(
θ
i
,ξ
i
)
}
N
t
i
=1
of model parameter to model output. Harnessing the input-output pairs
192
as training points, Gaussian process regression is used (Rasmussen, 2003) to create an emu-
193
lator, composed of a mean function and covariance function pair, where
G
GP
(
θ
)
≈ G
∞
(
θ
)
194
and
Σ
GP
≈
Σ
. Since the input-output pairs are subject to different realizations of the chaotic
195
system, the emulator mean approximates
G
∞
(
θ
)
rather than
G
(
θ
,ξ
)
(Cleary et al., 2021; Dun-
196
bar et al., 2020).
197
The variables of interest in the synthetic data
Y
∈
R
N
×
N
d
are correlated. The output
variables from the GCM forward model
G
(
θ
,ξ
)
are similarly correlated. The correlation be-
tween the GCM statistics results in nonzero off-diagonal covariance matrix elements. In or-
der to maintain a diagonal covariance in the GP emulator
Σ
GP
, we transform the GCM statis-
–7–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
tics into a decorrelated space using principal component analysis on
Σ
∈
R
N
×
N
(Cleary et
al., 2021). That is, we decompose the covariance matrix using the singular value decompo-
sition
Σ =
V D
2
V
ᵀ
,
(13)
with a matrix of principal component vectors (or singular vectors)
V
and the diagonal matrix
198
D
containing the square roots of the singular values
σ
i
. Often, in practical data assimilation
199
problems,
N > N
d
(Houtekamer & Zhang, 2016). In this case, the covariance matrix is rank
200
deficient, with
rank(Σ)
≤
min(
N,N
d
)
, with singular values
σ
2
i
= 0
for
i >
rank(Σ)
. Meth-
201
ods to decorrelate the data and model output with rank deficient covariance matrices are dis-
202
cussed in Appendix A.
203
The GP is trained using input-output pairs in the decorrelated space
{
θ
i
,
̃
G
(
θ
i
,ξ
i
)
}
N
t
i
=1
with the decorrelated space denoted by tildes,
̃
(
·
)
. The resulting input-output mapping is ap-
proximated as
̃
G
∞
(
θ
)
≈
N
(
̃
G
GP
(
θ
)
,
̃
Σ
GP
(
θ
))
.
(14)
The GP kernels used are a white-noise kernel added to an Automatic Relevance Determina-
204
tion (ARD) radial basis function kernel; further details provided in Dunbar et al. (2020); Cleary
205
et al. (2021). The GP hyperparameters are trained using the input-output pairs. The training
206
process results in a GP regression function
̃
G
GP
(
θ
)
which takes
θ
as input and emulates
̃
G
∞
(
θ
)
207
in a computationally efficient fashion.
208
2.4 Sample: Posterior sampling using MCMC
209
The Bayesian posterior distribution is approximated through MCMC sampling with the
trained GP emulator
̃
G
GP
(
θ
)
. The data
y
is normalized and transformed into the decorrelated
space
̃
y
, as in section 2.3. The Bayesian posterior distribution is approximated as (Stuart, 2010)
P
(
θ
|
̃
y
)
∝
P
(
̃
y
|
θ
)
P
(
θ
)
(15)
∝
exp
(
−
1
2
‖
̃
y
−
̃
G
GP
(
θ
)
‖
2
̃
Σ
GP
(
θ
)+
̃
∆
−
1
2
log det
(
̃
Σ
GP
(
θ
) +
̃
∆
)
)
P
(
θ
)
.
(16)
With a normal prior distribution governed by mean
θ
and covariance
Γ
θ
the resulting MCMC
objective function is
Φ
MCMC
(
θ
,
̃
y
) = exp
(
−
1
2
‖
̃
y
−
̃
G
GP
(
θ
)
‖
2
̃
Σ
GP
(
θ
)+
̃
∆
−
1
2
log det
(
̃
Σ
GP
(
θ
) +
̃
∆
)
−
1
2
‖
θ
−
θ
‖
2
Γ
θ
)
.
(17)
We use a random walk Metropolis MCMC algorithm. The number of MCMC samples is set
210
to
N
MCMC
= 200
,
000
with a burn-in of
10
,
000
.
211
3 Seasonal GCM uncertainty quantification
212
We perform numerical experiments with an idealized GCM with a seasonal cycle. The
213
GCM simulation setup is presented in section 3.1. Various climate statistics that we extract
214
from the GCM and use in the numerical experiments are discussed in section 3.2 and shown
215
in section 3.3.
216
3.1 Seasonal simulation setup
217
The GCM used in this study is based on the Geophysical Fluid Dynamics Laboratory’s
218
Flexible Modeling System (Frierson et al., 2006). The GCM simulates an idealized aquaplanet
219
with a homogeneous mixed-layer slab ocean bottom boundary condition with a depth of
1 m
.
220
The GCM has been used previously for simulations of the hydrological cycle over a range of
221
–8–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
climates (O’Gorman & Schneider, 2008) and to characterize seasonal variability in the trop-
222
ics (Bordoni & Schneider, 2008; Merlis et al., 2013a, 2013b; Kaspi & Schneider, 2011; Bischoff
223
& Schneider, 2014, 2016; Wei & Bordoni, 2018). The GCM is axisymmetric and statistically
224
cyclostationary. The spectral transform method is used in the horizontal directions, and finite
225
differencing in sigma coordinates is used in the vertical direction. The horizontal resolution
226
used is T21 with
N
φ
= 32
discrete latitude points on the transform grid, and
20
sigma lev-
227
els (
σ
p
=
p/p
s
, where
p
is the pressure and
p
s
is the local surface pressure).
228
A two-stream gray radiation scheme is used. The top-of-atmosphere (TOA) insolation
229
is prescribed and varies according to a seasonal cycle (Bordoni & Schneider, 2008; Wei & Bor-
230
doni, 2018). The diurnal cycle insolation variations are neglected, with a daily average inso-
231
lation applied. The longwave and shortwave optical thicknesses depend on the latitude and pres-
232
sure. The radiative effects of variations of atmospheric water vapor or clouds are neglected,
233
and therefore, water vapor and cloud feedbacks are not included in the GCM.
234
The convection is parameterized using a simplified quasi-equilibrium Betts-Miller (SBM)
scheme (A. K. Betts, 1986; A. Betts & Miller, 1986; Frierson, 2007). Vertical profiles of tem-
perature and humidity are used to calculate precipitation and associated temperature and hu-
midity changes through a relaxation to moist adiabatic reference profiles (Frierson, 2007). The
relaxation is included as a forcing to the temperature and humidity balances
∂T
∂t
+
···
=
S
T
−
f
T
T
−
T
ref
τ
(18)
∂q
∂t
+
···
=
S
q
−
f
q
f
T
q
−
q
ref
τ
,
(19)
where the dots represent the dynamical terms in the equations, and
S
T
and
S
q
represent the
235
temperature and specific humidity forcings aside from the convection scheme. The term
f
T
236
governs the spatiotemporal activation of the convection scheme, depends on the thermodynamic
237
state, and is a function of
z
. The term
f
q
modifies the specific humidity relaxation (O’Gorman
238
& Schneider, 2008). The reference temperature
T
ref
is a moist adiabat, chosen so that the con-
239
vection scheme conserves enthalpy integrated over vertical columns (Frierson, 2007; O’Gorman
240
& Schneider, 2008). The reference specific humidity
q
ref
is that which corresponds to a pre-
241
scribed relative humidity with respect to the moist adiabat
T
ref
. For the UQ experiments, our
242
focus are two parameters: the prescribed reference relative humidity parameter (
θ
RH
) and the
243
relaxation timescale parameter (
θ
τ
).
244
The GCM simulation starts from vernal equinox, and the year length is
T
yr
= 360 d
.
245
All climate statistics used for UQ are zonally averaged. We use the Betts-Miller convection
246
scheme with standard reference values of the parameters of
θ
†
= (
θ
†
RH
,θ
†
τ
) = (0
.
7
,
7200 s)
247
(Frierson, 2007; O’Gorman & Schneider, 2008). UQ of
θ
relies on a prior knowledge about
248
the convective parameters. We use wide prior distributions of the parameters, which enforce
249
physical constraints, such as
0
< θ
RH
≤
1
and
θ
τ
>
0
, but are otherwise uninformative
250
(Dunbar et al., 2020). The selected priors are
θ
RH
∼
Logit(
N
(0
,
1))
and
θ
τ
∼
Log(
N
(12h
,
(12h)
2
))
,
251
that is, normal distributions of logit- and log-transformed parameters. The parameter priors
252
are independent, although joint prior distributions could be used in future work.
253
3.2 Climate statistics used for UQ
254
We calibrate and perform UQ on the parameters of the convection scheme using more
255
and less informative time-averaged data from GCM simulations. The informative statistics, such
256
as the mid-tropospheric relative humidity, are chosen because they are strongly affected by the
257
choice of the convective parameters (e.g., O’Gorman and Schneider (2008); Schneider and O’Gorman
258
(2008); O’Gorman et al. (2011)). For comparison, we also choose less informative statistics
259
that are affected less by the choice of convective parameters, such as surface wind speeds. Nu-
260
merical UQ experiments using the differing degrees of information in the climate statistics will
261
illustrate the impact of incorporating higher frequencies in the climate statistics used for pa-
262
rameter estimation.
263
–9–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
Informative
Less informative
Variables
1. Relative humidity (
σ
p
= 0
.
5
)
2. Precipitation rate (mm/day)
3. Probability of
90
th
percentile precipitation
1. Precipitation rate (mm/day)
2. Surface wind speed (m/s)
Latitudes,
N
φ
32
32
N
y
96
64
N, T
= 360 d
(annual)
96
64
N, T
= 90 d
(seasonal)
384
256
Table 1: Summary of the climate data used. Calibration and UQ is performed using data that are
more or less informative about the convective scheme. The total size of the data used for UQ is
N
.
For the informative climate statistics, we choose three variables: the mid-tropospheric
264
(
σ
p
= 0
.
5
) relative humidity, the total precipitation rate, and, as a measure of precipitation
265
intensity, the probability that a daily precipitation total exceeds the latitude-dependent 90th pre-
266
cipitation percentile threshold from a long control simulation with the true parameters (Dunbar
267
et al., 2020). Since intense precipitation events are influenced by the convection scheme (O’Gorman
268
& Schneider, 2009a, 2009b), exceedances of precipitation over a high threshold are anticipated
269
to be informative about the convection scheme parameters. The three statistics, are evaluated
270
at each of the
N
φ
= 32
discrete latitudes, giving
96
total quantities of interest. For less in-
271
formative statistics, we use the zonally averaged precipitation rate and surface wind speed. As
272
with the informative statistics,
N
φ
= 32
discrete latitudes are considered, giving
64
total quan-
273
tities of interest (see Table 1 for a summary).
274
We nondimensionalize the GCM statistics (Eq. (A1)) with the median
y
c
, taken over lat-
275
itude and time, of each specified data type. There is one characteristic value for each data type,
276
e.g., one characteristic value for precipitation and a separate characteristic value for relative
277
humidity.
278
3.3 Synthetic data used for UQ
279
As synthetic data, we used GCM output with different initial conditions, and perturbed
280
with measurement error with covariance
∆
. The synthetic states are constructed by running
281
the seasonal GCM for
150
years at the true parameters
θ
†
. For each case of
T
= 90 d
and
282
T
= 360 d
, the GCM states are aggregated according to the method introduced in section
283
2.1. The informative synthetic statistics for
T
= 360 d
days are shown in Figure 2a–c, which
284
compares the
N
d
= 150
year ensemble average of the statistics with a
1
year sample sub-
285
ject to internal variability and perturbed with measurement noise
N
(0
,
∆)
. The synthetic data
286
for
T
= 90 d
are shown in Figure 2d–f. For the data averaged over
90 d
, a distinct seasonal
287
cycle emerges with, for example, relatively high precipitation and relative humidity in the north-
288
ern hemisphere summer (days 90–180). The less informative statistics of precipitation rate and
289
surface wind speed are shown for
T
= 360 d
and
T
= 90 d
in Figure 3.
290
4 GCM calibration and UQ results
291
4.1 Informative statistics
292
The convective parameters
θ
RH
and
θ
τ
are calibrated using EKI with the informative statis-
293
tics with
T
= 90 d
and
T
= 360 d
(see Table 1). The EKI calibration is performed with a
294
synthetic data sample
y
(Eq. 7). Since the synthetic data is subject to internal variability, the
295
corresponding calibration is influenced by the synthetic sample values. The synthetic data sam-
296
–10–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
(a)
(a
a
aaa
(a)
(b)
(c)
(d)
(e)
(f)
Figure 2: Informative synthetic data with (a–c)
T
= 360 d
and (d–f)
T
= 90 d
. The light solid
lines correspond to
y
∞
. The shaded regions correspond to a
95%
confidence interval around
y
∞
with covariance
Σ
+
∆
. The dark solid lines with circle markers correspond to a randomly se-
lected sample of the synthetic data
y
that has been subjected to measurement noise
N
(0
,
∆)
. Day
0
corresponds to vernal equinox and the Northern hemisphere seasons are provided.
(a)
(b)
(c)
(d)
Figure 3: Less informative synthetic data with (a, b)
T
= 360 d
and (c, d)
T
= 90 d
. The light
solid lines correspond to
y
∞
. The shaded regions correspond to a
95%
confidence interval around
y
∞
with covariance
Σ
+
∆
. The dark solid lines with circle markers correspond to a randomly
selected sample of the synthetic data
y
that has been subjected to measurement noise
N
(0
,
∆)
.
Day
0
corresponds to vernal equinox and the Northern hemisphere seasons are provided.
–11–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
(a)
(b)
Figure 4: Ensemble Kalman inversion convective parameter estimates using informative statistics.
Mode of (a) the estimated relative humidity convective parameter
θ
RH
compared to truth value
θ
†
RH
and of (b) the estimated relaxation time scale convective parameter
θ
τ
compared to truth
value
θ
†
τ
. The mean and standard deviation over independent instantiations with differing synthetic
samples are shown by the lines and shaded regions, respectively.
ple is randomly selected from the
150
years of historical data constructed with the true con-
297
vective parameters
θ
†
(see section 3.2). Calibration is performed with
10
independent synthetic
298
samples, and the results in this section are presented as means and standard deviations over
299
the
10
independent realizations. The cases with
T
= 90 d
and
T
= 360 d
are run with the
300
same
10
synthetic samples.
301
The modes of the EKI ensemble distributions are shown in Figure 4. The modes demon-
302
strate a reduction in estimation error of the true convective parameters when the data are ag-
303
gregated seasonally with
T
= 90 d
. The calibrated convective parameters using annually av-
304
eraged GCM statistics have higher mean bias and variability between the instantiations with
305
different synthetic data samples. For both integration timescales, the error associated with the
306
θ
τ
estimation was larger than the estimation error for
θ
RH
. The percent error of the average
307
mode at the last EKI iteration for
θ
RH
is
0
.
3%
and
0
.
8%
for
T
= 90 d
and
T
= 360 d
,
308
respectively. For
θ
τ
, the percent errors are
1
.
2%
and
11
.
5%
for
T
= 90 d
and
T
= 360 d
,
309
respectively. This indicates that the integration of the GCM states over an annual cycle filters
310
information from the resulting GCM statistics that is informative, especially about the relax-
311
ation timescale.
312
The normalized mean square error (MSE) for
θ
RH
at EKI iteration
(
n
)
is
ε
(
n
)
RH
=
1
θ
†
RH
√
√
√
√
1
N
ens
N
ens
∑
i
=1
(
θ
RH
,i
−
θ
†
RH
)
2
.
(20)
The error is computed similarly for
θ
τ
. For both the relative humidity and relaxation timescale
313
parameters,
T
= 90 d
reduces the parameter estimation MSE relative to
T
= 360 d
by about
314
a factor 2–3: The MSE for
θ
RH
is
0
.
006
and
0
.
015
for
T
= 90 d
and
T
= 360 d
, respec-
315
tively; for
θ
τ
, the MSE is
0
.
06
and
0
.
23
for
T
= 90 d
and
T
= 360 d
, respectively. The
316
seasonally aggregated data also have a smaller spread among the realizations, as visualized
317
by the uncertainty bands around the mean values. The reductions in the MSE standard devi-
318
ations were
43%
and
66%
for
θ
RH
and
θ
τ
, respectively.
319
The input-output pairs generated during EKI parameter calibration are used for training
320
the GP emulator. These pairs of parameters and model evaluations are taken from four 100-
321
member EKI iterations as well the initial 100-member ensemble drawn from the prior, giv-
322
ing a total of 500 input-output pairs for GP training. The results of this study were not sig-
323
nificantly different with more training points. The GCM outputs are mapped into a normal-
324
ized and decorrelated space according to the SVD of the internal variability covariance ma-
325
–12–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
(a)
(b)
(c)
(d)
(e)
(f)
Figure 5: Convective parameter posterior distributions computed using MCMC using a GP em-
ulator that was trained using data aggregated with an integration timescale of
T
=
360 d
and
T
= 90 d
. The contours correspond to
50%
,
75%
, and
99%
of the posterior distribution, the star is
the truth, and the circle is the average of the ensemble members after the last EKI iteration. Poste-
riors are shown for (a)
T
= 360 d
, decorrelated with full SVD, and (b) a 95%-energy truncation of
the SVD. (c) Posterior for
T
= 90 d
, decorrelated with 95%-energy truncated SVD. Panels (d–f)
same as (a–c) with the posterior distributions shown in the physical parameter space.
–13–
ESSOAr | https://doi.org/10.1002/essoar.10507655.1 | Non-exclusive | First posted online: Wed, 4 Aug 2021 10:53:11 | This content has not been peer reviewed.