{'fraction_non_alphanumeric': 0.03334964282219395, 'fraction_numerical': 0.015852823172521773, 'mean_word_length': 4.810986011600137, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Parameter uncertainty quantification in an idealized GCM with a seasonal cycle

Climate models are generally calibrated manually by comparing selected climate statistics, such as the global top-of-atmosphere energy balance, to observations. The manual tuning only targets a limited subset of observational data and parameters. Bayesian calibration can estimate climate model parameters and their uncertainty using a larger fraction of the available data and automatically exploring the parameter space more broadly. In Bayesian learning, it is natural to exploit the seasonal cycle, which has large amplitude, compared with anthropogenic climate change, in many climate statistics. In this study, we develop methods for the calibration and uncertainty quantification (UQ) of model parameters exploiting the seasonal cycle, and we demonstrate a proof-of-concept with an idealized general circulation model (GCM). Uncertainty quantification is performed using the calibrate-emulate-sample approach, which combines stochastic optimization and machine learning emulation to speed up Bayesian learning. The methods are demonstrated in a perfect-model setting through the calibration and UQ of a convective parameterization in an idealized GCM with a seasonal cycle. Calibration and UQ based on seasonally averaged climate statistics, compared to annually averaged, reduces the calibration error by up to an order of magnitude and narrows the spread of posterior distributions by factors between two and five, depending on the variables used for UQ. The reduction in the size of the parameter posterior distributions leads to a reduction in the uncertainty of climate model predictions.


Introduction
The objective of quantifying uncertainty in computational models arises in a wide range of applications, including weather and climate modeling [1], fluid dynamics [2], and energy systems [3]. Often, uncertainty associated with predictions from computational models is the result of processes that cannot be resolved on the computational grid, either due to computational complexity limitations [4] or due to uncertainty associated with the process itself [5]. In general circulation models (GCMs), primary uncertainties arise from the representation of subgrid-scale turbulence, convection, and cloud physics, which have a significant impact on the evolution of climate under rising greenhouse gases [6,7,8,9,10]. While clouds are associated with turbulence and convective updrafts with scales of O(10 m), modern climate simulations have a typical horizontal resolution of O(10 km)-O(100 km) [11,12]. Climate simulations rely on physically motivated parameterizations that model the effects of subgrid-scale processes such as clouds and turbulence on the resolved scales [13]. Such parameterizations come with parametric and structural uncertainties; quantifying these uncertainties and how they percolate into climate projections remains an outstanding challenge [5]. arXiv:2108.00827v1 [physics.ao-ph] 22 Jul 2021 Physical parameterizations have historically been individually developed and calibrated using data from isolated experiments [14,15,16,17,18]. They are further adjusted so that global models that incorporate them satisfy selected large-scale observational or physical constraints, such as a closed top-of-atmosphere energy balance or reproduction of the 20th-century global-mean temperature evolution [19,18,20]. Model calibration is usually done manually, focusing on a small subset of model parameters and exploiting only a fraction of the available observational data.
As a step toward automating and augmenting this process, here we further develop algorithms for model calibration and uncertainty quantification (UQ) that in principle allow models to learn from large datasets and that scale to highdimensional parameter spaces. In previous work, these algorithms have been demonstrated in simple conceptual models [21] and in a statistically stationary idealized GCM [22]. We take the next step and demonstrate how these algorithms can exploit seasonal variations, which for many climate statistics are large relative to the climate changes expected in the coming decades and contain exploitable information about the response of the climate system to perturbations [23,24].
Whereas numerical weather prediction assimilates spatiotemporally evolving trajectories of atmospheric states as initial conditions for forecasts (e.g., [25,26]), in climate modeling it is preferable to assimilate time-averaged climate statistics. This focuses the learning problem on quantities of interest in climate predictions (i.e., climate statistics, including higher-order statistics such as precipitation extremes), and it avoids the need to estimate uncertain atmospheric initial conditions on which trajectories of states depend [21]. Calibration and UQ of climate models on the basis of time-averaged statistics smooths the prediction-error based objective function and enables the use of data that have different resolution than the climate simulations [1,22]. Given the dual desires to avoid having to estimate atmospheric initial conditions, which are forgotten over about 2 weeks [27], and to exploit seasonal variations, it becomes natural to choose averaging timescales between around 30 and 90 days. Such averaging timescales are the focus of this study.
We consider the calibration and UQ of convective parameters in an idealized GCM with seasonally varying insolation [28,29,30]. We develop an extension of the calibrate-emulate-sample (CES) Bayesian learning methodology [21] to enable the use of statistics computed from a non-stationary statistical state. The qualitative and quantitative impacts of the time-averaging length of the climate statistics on the parameter calibration and UQ are assessed. GCMs are generally tuned in situations which have low parameter identifiability based on available climate data [18]. We perform numerical experiments with observable climate statistics that are informative about the convective parameters we wish to calibrate, but are not in any simple and direct way related to them. We also explore less informative climate statistics, which highlight the benefits of incorporating seasonal variations in the climate statistics that are being exploited for UQ.
The remainder of this paper is organized as follows. In section 2, the Bayesian learning methods for time-dependent problems are introduced. The numerical details of the seasonally forced GCM and UQ experiments are introduced in section 3. The calibration and UQ results are shown in section 4, and conclusions are provided in section 5.

Uncertainty quantification methods
The goal of this study is to estimate the probability distributions associated with model parameters θ that are used by a GCM and about which only imprecise prior information is known. We will consider a Bayesian approach to the estimation of the probability distributions of model parameters θ, where we seek P(θ|y), the conditional probability distribution of parameters θ given observed data y. The GCM is a computationally expensive numerical model that evolves climate states in time. From the GCM, we extract statistical information that is denoted by G(θ). Here, G(θ) includes the numerical integration of the GCM states and the aggregation of relevant climate statistics in time. Although the data y will, in general, not provide direct information about θ, we can estimate P(θ|y) using Bayesian learning by comparing the GCM outputs G(θ) with data y.
Standard methods for Bayesian learning include Markov Chain Monte Carlo (MCMC) [31], which typically requires O(10 5 ) forward model evaluations of G(θ) to sample the posterior distribution [32]. Instead, we perform calibration and UQ using the recently developed CES methodology [21], which consists of three steps: (1) Ensemble Kalman processes [33,34] are used to calibrate parameters and to generate input-output pairs of the mapping θ → G(θ); (2) Gaussian process (GP) regression is used to train an emulator G GP (θ) of the mapping θ → G(θ) using the training points generated in the calibration step; and (3) MCMC sampling with the computationally efficient GP emulator G GP (θ) rather than the expensive forward model G(θ) is used to estimate the posterior distribution P(θ|y). The CES methodology has previously been used for calibration and UQ of parameters in simple model problems such as Darcy flow and Lorenz systems [21] and for convective parameters in a statistically stationary GCM [22]. More recent methodological developments by [35] enabled the CES framework to perform simultaneous UQ on O(1000) parameters using deep neural network-based emulation and MCMC suited to high-dimensional spaces [36,37]. A schematic of the CES methodology is shown in Figure 1.

EKI/EKS using
Inputs : Figure 1: Schematic of the calibrate-emulate-sample (CES) methodology to estimate model parameters θ. 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, uncorrelated 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 posterior distribution P(θ|ỹ). The objective functions for calibration and sampling are denoted by Φ(θ, y) and Φ MCMC (θ,ỹ), respectively.
The Bayesian learning methodology used in this study is described in the following sections. In section 2.1, the inverse problem of estimating P(θ|y) in a setting with a periodic cycle is introduced. In section 2.2, the ensemble Kalman process calibration method is outlined. Section 2.3 introduces the GP emulation in an uncorrelated transformed space, obtained by a singular value decomposition (principal component analysis) of the noise covariance matrix. Section 2.4 describes how the Bayesian posterior distribution is approximated using the GP emulator.

Seasonal GCM inverse problem
With fixed insolation [28,29], GCM statistics 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 cyclostationary, with a dependence on the ordinal day. The data we use are constructed accounting for the seasonally varying boundary conditions, as the time-averaged GCM statistics Here, H(t; θ, ξ) ∈ R Ny represents the output states of the GCM, depending on time t. The number of states being measured is N y . The integration length is specified as T . The initial conditions for the GCM integration are represented by ξ, which depend on the initial time t. The time-averaged statistics are G T (θ, ξ, t) ∈ R Ny .
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 set of time-averaged statistics where j = 1 through N s is an index representing the time of year in the GCM simulation, starting from j = 1 at vernal equinox (time t 0 ). The integration windows are non-overlapping. The length of the year in the GCM is T yr = 360 d, and the resulting number of non-overlapping time-averaged statistics are N s = T yr /T . In this framework, the number of batches of statistics N s is a design parameter set by the integration timescale T . For T = 90 d, the data are aggregated seasonally, with N s = 4. For T = 360 d, the statistics are averaged over the full year in the GCM, corresponding to annually averaged climate statistics, and N s = 1.
The size of y and G is R N , where N = N s · N y . The resulting inverse problem relating parameters and data is where η ∼ N (0, ∆) is a realization of normal measurement error with zero mean and covariance matrix ∆. We generate synthetic data and forward model data starting from the same initial conditions ξ in Eq. (5). Since the states of the GCM depend on the ordinal day, the averaging operation in Eq. (2) must be consistently aligned in ordinal day between the synthetic data and GCM outputs.
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 different initial conditions ξ i is In the limit of infinite realizations of the climate statistics, the dependence of the ensemble average G ∞ (θ) on the initial conditions trends to zero by the central limit theorem. The central limit theorem applies since G(θ, ξ i ) contains a full year of data by construction (Eq. 3) and is therefore a statistically stationary object with respect to the cycle. But due to computational limitations, only a finite number of ensemble realizations is available in practice. For the synthetic data, we can similarly define y ∞ ; however, y ∞ is generally not accessible from observations.
The atmospheric initial condition is of minimal practical value in the climate model setting, 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 estimation 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, where Σ is the internal variability covariance matrix of the GCM [21]. The realizations G(θ, ξ) and y are both subject to the internal variability 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 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 dimension 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 where · A = A −1/2 · 2 . The likelihood is [38] P(θ|y) ∝ exp(−Φ(θ, y)).
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 error. To emulate measurement error, we add Gaussian noise to the GCM output, with zero mean and covariance matrix Here, the noise standard deviation is defined such that where C m = 0.1 caps the maximum measurement error standard deviation to 10% of the mean data values and C = 0.2 controls how close the noise-perturbed data can come within physical boundaries ∂Ω i of the data (e.g., to keep relative humidities between 0 and 1 with high probability).

Calibrate: Ensemble Kalman Inversion
The first stage of the CES methodology is to calibrate the parameters θ of the model based on data y. We perform calibration with independent realizations of G(θ i , ξ i ), viewed as noisy approximations of G ∞ (θ i ). Calibration is performed using ensemble Kalman methods, which demonstrate theoretical success, in idealized problems, and empirical success, in complex problems, to optimize parameters under such noise [39].
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 [40]. The EnKF has been used widely for derivative-free state estimation in numerical weather prediction (e.g., [41]) and model-based control (e.g., [42]). Ensemble Kalman methods for Bayesian inversion were introduced by [43] and [44]. These methods provably 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) [45] algorithm was introduced for classical, optimization-based inversion. EKI generally drives the ensemble members toward consensus near the optimal solution of the inverse problem [33]. The parameter update of ensemble member m at iteration step n is where C θG is the empirical cross-covariance between the parameters and the model outputs, and C GG is the empirical covariance of the model outputs.
EKI is guaranteed to find the optimal parameters in linear problems [33] and has empirical success in nonlinear problems [22]. Several other approaches exist for estimation in nonlinear problems. [46] and [47] developed iterative EnKF methods which improve state estimation in strongly nonlinear problems. However, the distribution of the ensemble does not converge to the posterior distribution in the limit of infinite members for nonlinear problems [48,49,50,34], necessitating the emulation and sampling described in the following sections.
The number of ensemble members and EKI iterations are hyperparameters; they are set to standard values of N ens = 100 and N it = 5 in this study [33,22,21]. Each ensemble member is run through the GCM, initialized from the same initial conditions. A spin-up period of one year (360 days) is run before the statistics are computed using Eq. (2), to ensure the forward model realizations are subject to differing instantiations of internal variability in the chaotic GCM system.

Emulate: Gaussian process emulation
We emulate the mapping from parameters to model output using a machine learning method that enables the rapid execution of the mapping, compared to the computationally expensive forward model. The calibration stage (section of model parameter to model output. Harnessing the input-output pairs as training points, Gaussian process regression is used [51] to create an emulator, composed of a mean function and covariance function pair, where G GP (θ) ≈ G ∞ (θ) and Σ GP ≈ Σ. Since the input-output pairs are subject to different realizations of the chaotic system, the emulator mean approximates G ∞ (θ) rather than G(θ, ξ) [21,22].
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 between the GCM statistics results in nonzero off-diagonal covariance matrix elements. In order to maintain a diagonal covariance in the GP emulator Σ GP , we transform the GCM statistics into a decorrelated space using principal component analysis on Σ ∈ R N ×N [21]. That is, we decompose the covariance matrix using the singular value decomposition with a matrix of principal component vectors (or singular vectors) V and the diagonal matrix D containing the square roots of the singular values σ i . Often, in practical data assimilation problems, N > N d [41]. In this case, the covariance matrix is rank deficient, with rank(Σ) ≤ min(N, N d ), with singular values σ 2 i = 0 for i > rank(Σ). Methods to decorrelate the data and model output with rank deficient covariance matrices are discussed in A.
The GP is trained using input-output pairs in the decorrelated space with the decorrelated space denoted by tildes,(·). The resulting input-output mapping is approximated as The GP kernels used are a white-noise kernel added to an Automatic Relevance Determination (ARD) radial basis function kernel; further details provided in [22,21]. The GP hyperparameters are trained using the input-output pairs. The training process results in a GP regression functionG GP (θ) which takes θ as input and emulatesG ∞ (θ) in a computationally efficient fashion.

Sample: Posterior sampling using MCMC
The Bayesian posterior distribution is approximated through MCMC sampling with the trained GP emulatorG GP (θ). The data y is normalized and transformed into the decorrelated spaceỹ, as in section 2.3. The Bayesian posterior distribution is approximated as [52] P(θ|ỹ) ∝ P(ỹ|θ)P(θ) With a normal prior distribution governed by mean θ and covariance Γ θ the resulting MCMC objective function is We use a random walk Metropolis MCMC algorithm. The number of MCMC samples is set to N MCMC = 200, 000 with a burn-in of 10, 000.

Seasonal GCM uncertainty quantification
We perform numerical experiments with an idealized GCM with a seasonal cycle. The GCM simulation setup is presented in section 3.1. Various climate statistics that we extract from the GCM and use in the numerical experiments are discussed in section 3.2 and shown in section 3.3.

Seasonal simulation setup
The GCM used in this study is based on the Geophysical Fluid Dynamics Laboratory's Flexible Modeling System [28]. The GCM simulates an idealized aquaplanet with a homogeneous mixed-layer slab ocean bottom boundary condition with a depth of 1 m. The GCM has been used previously for simulations of the hydrological cycle over a range of climates [29] and to characterize seasonal variability in the tropics [30,53,54,55,56,57,58]. The GCM is axisymmetric and statistically cyclostationary. The spectral transform method is used in the horizontal directions, and finite differencing in sigma coordinates is used in the vertical direction. The horizontal resolution used is T21 with N φ = 32 discrete latitude points on the transform grid, and 20 sigma levels (σ p = p/p s , where p is the pressure and p s is the local surface pressure).
A two-stream gray radiation scheme is used. The top-of-atmosphere (TOA) insolation is prescribed and varies according to a seasonal cycle [30,58]. The diurnal cycle insolation variations are neglected, with a daily average insolation applied. The longwave and shortwave optical thicknesses depend on the latitude and pressure. The radiative effects of variations of atmospheric water vapor or clouds are neglected, and therefore, water vapor and cloud feedbacks are not included in the GCM.
The convection is parameterized using a simplified quasi-equilibrium Betts-Miller (SBM) scheme [59,60,61]. Vertical profiles of temperature and humidity are used to calculate precipitation and associated temperature and humidity changes through a relaxation to moist adiabatic reference profiles [61]. The relaxation is included as a forcing to the temperature and humidity balances where the dots represent the dynamical terms in the equations, and S T and S q represent the temperature and specific humidity forcings aside from the convection scheme. The term f T governs the spatiotemporal activation of the convection scheme, depends on the thermodynamic state, and is a function of z. The term f q modifies the specific humidity relaxation [29]. The reference temperature T ref is a moist adiabat, chosen so that the convection scheme conserves enthalpy integrated over vertical columns [61,29]. The reference specific humidity q ref is that which corresponds to a prescribed relative humidity with respect to the moist adiabat T ref . For the UQ experiments, our focus are two parameters: the prescribed reference relative humidity parameter (θ RH ) and the relaxation timescale parameter (θ τ ).
The GCM simulation starts from vernal equinox, and the year length is T yr = 360 d.  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 . [61,29]. UQ of θ relies on a prior knowledge about the convective parameters. We use wide prior distributions of the parameters, which enforce physical constraints, such as 0 < θ RH ≤ 1 and θ τ > 0, but are otherwise uninformative [22]. The selected priors are θ RH ∼ Logit (N (0, 1)) and θ τ ∼ Log (N (12h, (12h) 2 )), that is, normal distributions of logit-and log-transformed parameters. The parameter priors are independent, although joint prior distributions could be used in future work.

Climate statistics used for UQ
We calibrate and perform UQ on the parameters of the convection scheme using more and less informative time-averaged data from GCM simulations. The informative statistics, such as the mid-tropospheric relative humidity, are chosen because they are strongly affected by the choice of the convective parameters (e.g., [29,62,63]). For comparison, we also choose less informative statistics that are affected less by the choice of convective parameters, such as surface wind speeds. Numerical UQ experiments using the differing degrees of information in the climate statistics will illustrate the impact of incorporating higher frequencies in the climate statistics used for parameter estimation.
For the informative climate statistics, we choose three variables: the mid-tropospheric (σ p = 0.5) relative humidity, the total precipitation rate, and, as a measure of precipitation intensity, the probability that a daily precipitation total exceeds the latitude-dependent 90th precipitation percentile threshold from a long control simulation with the true parameters [22]. Since intense precipitation events are influenced by the convection scheme [64,65], exceedances of precipitation over a high threshold are anticipated to be informative about the convection scheme parameters. The three statistics, are evaluated at each of the N φ = 32 discrete latitudes, giving 96 total quantities of interest. For less informative statistics, we use the zonally averaged precipitation rate and surface wind speed. As with the informative statistics, N φ = 32 discrete latitudes are considered, giving 64 total quantities of interest (see Table 1 for a summary).
We nondimensionalize the GCM statistics (Eq. (21)) with the median y c , taken over latitude and time, of each specified data type. There is one characteristic value for each data type, e.g., one characteristic value for precipitation and a separate characteristic value for relative humidity.

Synthetic data used for UQ
As synthetic data, we used GCM output with different initial conditions, and perturbed with measurement error with covariance ∆. The synthetic states are constructed by running the seasonal GCM for 150 years at the true parameters

Informative statistics
The convective parameters θ RH and θ τ are calibrated using EKI with the informative statistics with T = 90 d and T = 360 d (see Table 1). The EKI calibration is performed with a synthetic data sample y (Eq. 7). Since the synthetic  : 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. data is subject to internal variability, the corresponding calibration is influenced by the synthetic sample values. The synthetic data sample is randomly selected from the 150 years of historical data constructed with the true convective parameters θ † (see section 3.2). Calibration is performed with 10 independent synthetic samples, and the results in this section are presented as means and standard deviations over the 10 independent realizations. The cases with T = 90 d and T = 360 d are run with the same 10 synthetic samples.
The modes of the EKI ensemble distributions are shown in Figure 4. The modes demonstrate a reduction in estimation error of the true convective parameters when the data are aggregated seasonally with T = 90 d. The calibrated convective parameters using annually averaged GCM statistics have higher mean bias and variability between the instantiations with different synthetic data samples. For both integration timescales, the error associated with the θ τ estimation was larger than the estimation error for θ RH . The percent error of the average mode at the last EKI iteration for θ RH is 0.3% and 0.8% for T = 90 d and T = 360 d, respectively. For θ τ , the percent errors are 1.2% and 11.5% for T = 90 d and T = 360 d, respectively. This indicates that the integration of the GCM states over an annual cycle filters information from the resulting GCM statistics that is informative, especially about the relaxation timescale.
The normalized mean square error (MSE) for θ RH at EKI iteration (n) is The error is computed similarly for θ τ . For both the relative humidity and relaxation timescale parameters, T = 90 d reduces the parameter estimation MSE relative to T = 360 d by about a factor 2-3: The MSE for θ RH is 0.006 and 0.015 for T = 90 d and T = 360 d, respectively; for θ τ , the MSE is 0.06 and 0.23 for T = 90 d and T = 360 d, respectively. The seasonally aggregated data also have a smaller spread among the realizations, as visualized by the uncertainty bands around the mean values. The reductions in the MSE standard deviations were 43% and 66% for θ RH and θ τ , respectively.
The input-output pairs generated during EKI parameter calibration are used for training the GP emulator. These pairs of parameters and model evaluations are taken from four 100-member EKI iterations as well the initial 100-member ensemble drawn from the prior, giving a total of 500 input-output pairs for GP training. The results of this study were not significantly different with more training points. The GCM outputs are mapped into a normalized and decorrelated space according to the SVD of the internal variability covariance matrix Σ, truncated to contain 95% of the energy in the singular values. For T = 90 d, the truncation contains k = 47 singular values, while T = 360 d retains k = 23 (see B). A scalar-valued GP is trained for each of the k outputs. For comparison, we also ran the emulation and sampling steps for T = 360 d without SVD truncation.
The posterior distributions are approximated using MCMC sampling with the trained GP emulators. Posterior distributions for the two convective parameters for a randomly selected synthetic data sample are shown in Figure 5. In Figure  5a, the posterior distribution was computed with normalized input data, aggregated annually with T = 360 d, but the SVD was not truncated before emulation and sampling. The posterior resulting from the same initial ensemble and synthetic data sample but with SVD truncation is shown in Figure 5b. The truncation of the SVD smooths and inflates 0.5% ± 0.1% 1.3% ± 0.3% Less informative data 5.5% ± 0.9% 26.0% ± 7.0% Table 2: The ratio (%) of the area occupied by the posterior to the area occupied by the prior. The area of the convex hull containing 75% of the posterior mass for each value of T and data type is normalized by the area containing 75% of the prior mass. The means and standard deviations over independent synthetic data realizations are provided. Ten and four realizations are used for the informative and less informative data, respectively. the posterior distribution. While the posteriors are qualitatively similar, they differ quantitatively. The true parameters are outside the region containing 99% of the posterior mass for the untruncated case but are within the region containing 75% of the posterior mass for the truncated case. Further investigation of the influence of SVD regularization on the posteriors is in B. Hereafter, all posterior results will focus on the posterior distributions estimated using emulation and sampling in the truncated SVD space (see A). The convective parameter posterior distribution using seasonally aggregated data with T = 90 d with the SVD truncated at 95% of the total energy is shown in Figure 5c. As in the calibration stage, the posterior mode for T = 90 d is closer to the true parameters than for T = 360 d. Compared to the posterior with T = 360 d, the T = 90 d posterior is also more compact.
We quantify the size of the posterior distributions for the two integration timescales by computing the two-dimensional area of the convex hull containing 75% of the posterior mass (middle contour level in Figure 5). The area containing 75% of the posterior mass for each value of T is normalized by the area containing 75% of the prior mass for reference. The percentages of the prior area taken up by the posteriors for T = 90 d and T = 360 d, averaged over the 10 CES instantiations, are shown in Table 2. The areas are computed in the logit(θ RH ) − log(θ τ /1 s) transformed space, in which the sampling is performed. The area of the posterior distribution resulting from the seasonal integration is (a) (b) Figure 6: Parameter calibration performed with EKI using less informative statistics. Mean square error of (a) the estimated relative humidity parameter θ RH compared to truth value θ † RH and of (b) the estimated relaxation time scale θ τ compared to truth value θ † τ . Solid lines are T = 360 d and dashed lines are T = 90 d.
approximately a factor 2 smaller than the posterior from the annually integrated data. However, compared to the wide prior, both posteriors are tight, with the posterior areas occupying around 1% of the prior area.

Less informative statistics
In many UQ applications, intuition about which quantities of interest will lead to improved parameter estimation is not available before experimentation. While the idealized GCM and convective parameterization used in this study have a rich set of investigations to provide prior knowledge about the relationship between θ and G(θ), we also consider UQ of the parameters using less informative GCM statistics. To do so, we use the precipitation rate and surface wind speed, zonally averaged and integrated with T = 360 d and T = 90 d (Figure 3). Four instantiations of CES with differing synthetic data samples are run for each integration length of T = 90 d and T = 360 d.
The MSE of the EKI parameter calibration with the less informative GCM statistics are shown in Figure 6 for T = 90 d and T = 360 d. With precipitation rate and surface wind speed as the basis for UQ, the annually averaged data result in significantly higher parameter estimation error than seasonally aggregated data. The T = 360 d MSEs are 0.26 and 1.56 for θ RH and θ τ , respectively. The T = 90 d MSEs are 0.02 and 0.28 for θ RH and θ τ , respectively. The ensemble means and modes of the two parameters are shown in Figure 7. In C, the calibration is performed by considering the seasonal data sequentially, rather than collectively, to restrict the data size for both T = 90 d and T = 360 d to N = 64.
The results demonstrate that the reductions in MSE are the result of the incorporation of seasonal information in the data, rather than the dimensionality of the data space.
Posterior distributions computed from CES with GP training data generated by EKI for T = 360 d and T = 90 d are shown in Figure 8. We use SVD truncation for both T = 360 d and T = 90 d (A). The posterior for T = 360 d is significantly larger than for T = 90 d, indicating that the annually averaged precipitation and surface wind speed do not provide substantive information regarding the convective parameters. As with EKI for T = 360 d, the posterior distribution mode has substantial error when compared to the true parameters. For T = 90 d, the uncertainty for θ RH estimated with the less informative data has collapsed to a similar order as the uncertainty estimated with informative data. While the uncertainty is larger for θ τ when estimated using the less informative data, the overall posterior area is reasonably collapsed, as it is of similar order to the T = 360 d posterior with informative data. These results suggest that seasonally averaged precipitation rate and surface wind speed are sufficient statistics to estimate convective parameters while annually averaged precipitation and wind speed introduce higher estimation error and parameter uncertainty.
As with the informative GCM statistics, we compare the sizes of the posterior distributions by computing the twodimensional area of the convex hull containing 75% of the posterior mass. The posterior areas, normalized by the prior area, averaged over the four CES instantiations, are shown in Table 2

Prediction experiments
To demonstrate the effect of parametric uncertainty in the GCM on climate predictions, we draw samples of 100 parameter pairs from the posterior distributions resulting from the less informative statistics of precipitation rate and surface wind speed (Figure 8). The samples of parameter pairs using seasonal integration (T = 90 d) and annual integration (T = 360 d) in the UQ are shown in Figure 9. With each sample of parameters, we produce climate quantities of interest, averaged over a 20 year period. We compute ensemble statistics over the outputs from the samples (prediction uncertainty) and compare these with a 20 year simulation with the true parameters θ † .
The relative humidity at σ p = 0.5, the precipitation rate, and the intense precipitation probability (the probability of exceeding a 90th percentile threshold averaged over the 20 year simulation at θ † ) are shown in Figure 10. The posterior distribution was estimated using (seasonally averaged) precipitation rate data directly, but neither relative humidity nor intense precipitation probability were used. The mean absolute percent errors for the mean prediction values, averaged across latitudes, for the T = 90 d posterior are 1.2%, 0.6%, and 1.4%, for the relative humidity, precipitation rate, and intense precipitation probability, respectively. For the predictions with the T = 360 d posterior, the mean absolute percent errors are 4.7%, 1.8%, and 21.5%.
The prediction uncertainty in each climate quantity of interest is indicated by a 95% confidence interval, estimated as the values of the 2.5 th and 97.5 th percentiles of the variables of interest at each latitude over the convective parameter pairs (Figure 9). For both T = 90 d and T = 360 d, there is limited prediction uncertainty in the mean precipitation rate. The widths of the confidence intervals are 6% and 4% of the mean precipitation rate, averaged across latitudes, for T = 360 d and T = 90 d, respectively. The prediction uncertainties in relative humidity and intense precipitation probability are significantly higher for UQ performed with annually averaged data (T = 360 d) compared to UQ with seasonally averaged data (T = 90 d). The T = 90 d posterior reduces the size of the 95% confidence interval for the relative humidity by 70%, averaged across latitudes, compared to T = 360 d. The reduction in the size of the 95% confidence interval for the intense precipitation probability is 9%.
We also performed idealized global-warming prediction experiments. To represent global warming in the idealized GCM, the longwave opacity in the atmosphere is increased by 50%, as in [29] which results in a global-mean surface air temperature increase from 284 K in the control climate to 292 K in the warm climate. We accumulated GCM statistics over a 20 year window after a spinup of one year, running simulations with the true GCM parameters θ † with increased longwave opacity for comparison. We use the same 100 parameter pairs drawn from the posterior distributions estimated from the less informative climate statistics.
The idealized global warming results are shown in Figure 11. relative humidity in a warmed climate by 75% compared to T = 360 d, averaged across latitudes. The reduction in the size of the 95% confidence interval for the intense precipitation probability in a warmed climate is 35%. These results demonstrate that the incorporation of the seasonal cycle in the UQ reduces the parametric uncertainty in climate predictions in this model.

Conclusions
We performed calibration and UQ of a convective parameterization in a seasonally forced idealized GCM. While GCMs are typically tuned using annually and globally averaged climate statistics [18], our results demonstrate, in an idealized setting, the qualitative, quantitative and systematic refinement of parameter distributions through the incorporation of seasonal information. Performing parameter calibration with seasonally averaged data significantly reduced the error associated with the estimated parameters compared to calibration using annually averaged data. The posterior distributions resulting from the Bayesian UQ with seasonally averaged data were reduced in size. One measure of success for Bayesian UQ, is to capture the true parameters within a high mass region of the posterior distribution (demonstrated in our results). We are also interested in the shape and size of the posterior distribution, as this provides valuable information about parameter correlations and uncertainty, with respect to observed data. We demonstrate that choosing suitable data for UQ which reduce the size of the posterior distribution can lead to significant reductions of climate prediction uncertainty.
The impact of incorporating additional frequency content is pronounced when the climate statistics used for UQ are less informative about the parameterizations. Such situations often occur in the climate modeling setting where a number of parameterizations may simultaneously and nonlinearly influence a quantity of interest and where it is not always clear which climate statistics should be used to calibrate an unknown GCM parameter.
To enable UQ in the climate model setting, we used the calibrate-emulate-sample (CES) methodology [21], which enables the efficient and accurate estimation of Bayesian posterior distributions of parameters from noisy climate data. CES uses gradient-free optimization to calibrate parameters and generate parameter-data pairs, Gaussian process regression to emulate the parameter-to-data mapping, and Markov Chain Monte Carlo (MCMC) to sample from the posterior distribution. The emulation and sampling is performed in a decorrelated dataspace through a transformation based on the principal component analysis of the noise covariance matrix. In this study, we modified the original CES methodology by first normalizing the data to ensure all statistics are weighted equally in the UQ and second by regularizing the covariance matrix to enable the UQ of applications with ill-conditioned or rank deficient covariance matrices. We quantified the impact of the normalization and regularization, with the regularization smoothing the posterior and slightly increasing its size.
Beyond convection, Earth system models rely on a number of parameterizations of cyclostationary processes, including models of carbon accumulation and storage [66] and of atmospheric boundary layer (ABL) turbulence which parameterize unclosed subgrid processes [67]. Carbon storage is inherently seasonal [68]. The ABL has distinct seasonal and diurnal variations [69,70]. We expect our results to be relevant to calibration and UQ for statistical variations on different timescales, beyond the seasonal cycle. For example, based on our findings presented here, we anticipate that incorporating the diurnal variation of ABL turbulence statistics may improve the UQ of subgrid scale turbulence models in GCMs. More broadly, we anticipate the findings of this study to be relevant to UQ in problems with multi-scale temporal dynamics.
The selection of the optimal aggregation timescale for climate statistics and choice of objective function remain open questions. In the context of parameter estimation, the integration timescale used to generate time-averaged statistics is a hyperparameter that dictates the trade-off between the frequency content and signal-to-noise ratio of the climate statistics. In the limiting case of T → 0, the influence of the initial conditions on the state is pronounced, as in weather forecasting [41]. Conversely, as T → ∞, all frequency content is removed from the statistics, which, as this study demonstrates, can adversely impact parameter estimation. We anticipate that the selection of the filter timescale T may be problem and parameterization specific, as it relates to the question of parameter identifiability. In this study, we selected seasonal averages because of the large amplitude of the seasonal cycle in many climate statistics and indications that seasonal variations are informative about the climate change response in many climate variables [24]. Future work should develop a more generalized approach for the selection of T based on the frequency content of the quantities of interest.

A Data and model output decorrelation for rank deficient problems
To facilitate the transformation into the uncorrelated space with rank deficient covariance matrices, the singular value decomposition/principal component analysis is truncated. Since the statistical quantities of interest have a range of magnitudes, the data are first normalized, and then the singular value decomposition is truncated as a form of regularization [71]. A detailed investigation of the influence of the normalization and truncation on the posteriors resulting from CES is performed in B.
The data used for UQ are provided concurrently to the CES pipeline in the concatenated vector y. In this framework, the statistical quantities in the data y may have a range of magnitudes, and normalization is required before regularization (e.g., [72]). The data are normalized by a characteristic value associated with each individual data type y c ∈ R N . The normalized data are where ( * ) denotes normalized data. The characteristic values in y c are described in the application of CES to the seasonally forced GCM in section 3. Each column of Y is normalized using Eq. (21) to yield the normalized data matrix Y * . The normalized covariance is Σ * = cov(Y * ). The SVD transformation is performed using covariance matrix Σ * and data y * .
The SVD is truncated in order to account for rank deficient or ill-conditioned covariance matrices [71]. The truncated SVD is defined for covariance Σ * .., σ * k ), truncated at the k ≤ rank(Σ * ) singular value. The truncated singular vector matrix is V k = [v 1 , ..., v k ] where v i is the singular vector corresponding to the singular value σ * 2 i . The truncated SVD space is given bỹ Parameter space Original Normalized Truncated Tikhonov Informative data 0.5% 0.5% 0.8% 0.9% Table 3: The ratio (%) of the area occupied by the posterior to the area occupied by the prior. The area of the convex hull containing 75% of the posterior mass for each regularization method is normalized by the area containing 75% of the prior mass.
In this study, the truncation location is selected as the lowest value of k such that k i=1 σ * 2 i ≥ 0.95 N i=1 σ * 2 i , so that 95% of the variance is retained. The GP output can be mapped to the original normalized space using Finally, the GP output can be transformed into the dimensional space through G GP = G * GP y c and Σ GP = Σ * GP y c y c , where denotes pointwise multiplication.
B Sensitivity of the posterior to covariance regularization In section 2.3, the SVD is truncated as a form of regularization. Here, we detail the sensitivity of the CES posterior distributions to the SVD regularization using informative GCM statistics (see Table 1) and annually averaged data with T = 360 d, for which the data covariance matrix is full rank. The data are normalized, as discussed in A. The cumulative sum of the singular values is shown in Figure 12. Truncation at 95% of the singular value energy corresponds to k = 23. Regularization can also be performed using Tikhonov regularization, which inflates the diagonal elements [71,73]. Following [71], regularization in the form of diagonal inflation Λ = diag(λ 2 ) is added to Σ where λ = (σ 3 k σ k+1 ) 1/4 . The SVD is performed such that  Table 3. The T = 360 d posterior distributions are shown in Figure 13. The regularization slightly increases the posterior area and smooths the posterior roughness.

C Parameter estimation with filtering and smoothing
It is often of practical interest to estimate parameters in problems that are not statistically stationary, such as the seasonally forced GCM with time-dependent boundary conditions. In problems which are not statistically stationary (e.g., with time-evolving boundary conditions), all requisite data for parameter UQ may not be available at the initialization of the estimation. This arises in the EnKF setting for state estimation (e.g., [41]).
We can instead pose EKI as a filtering approach for parameter estimation where the synthetic data are collected sequentially in a time-dependent problem. This differs from the definition of the data in section 2.1, where the data from each season are collected together into y (termed smoothing). For the filtering approach, the data is observed and averaged over time t j → t j + T in y T,j , where j indicates the time of the year in the GCM simulation. The EKI update is performed using y T,j and G T,j (θ (k) ) to update θ (k+1) , and this process is repeated for the specified number of EKI iterations. The smoothing approach leverages the data averaged in each season concurrently (y and G(θ (k) )) whereas the filtering approach uses the data averaged in each season sequentially (y T,j and G T,j (θ (k) )), as it becomes available.
The parameter calibration results comparing filtering and smoothing are shown in Figure 14. For T = 360 d, filtering and smoothing are identical. For T = 90 d, the filtering approach leverages the data from one season at a time to update the parameters θ. The final parameter errors for the T = 90 d filtering and smoothing approaches are similar, and are both significantly less than the error for the T = 360 d case. This result also demonstrates that the improvement in parameter estimates with T = 90 d, compared to T = 360 d, is a result of increased temporal information, rather than the dimensionality of the data space alone. Depending on the application, the filtering approach may be beneficial for computational efficiency, since it reduces the length of the forward model simulations, instead of requiring temporal integration over the timescale corresponding to an impactful low-frequency cycle.
(a) (b) Figure 14: Parameter calibration performed with EKI using less informative statistics. Mean square error 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 θ † τ . Solid lines are T = 360 d and dashed lines are T = 90 d.