{'fraction_non_alphanumeric': 0.03248203715423849, 'fraction_numerical': 0.00973530397230185, 'mean_word_length': 5.120186830713147, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improving seasonal forecast using probabilistic deep learning

The path toward realizing the potential of seasonal forecasting and its socioeconomic benefits depends heavily on improving general circulation model based dynamical forecasting systems. To improve dynamical seasonal forecast, it is crucial to set up forecast benchmarks, and clarify forecast limitations posed by model initialization errors, formulation deficiencies, and internal climate variability. With huge cost in generating large forecast ensembles, and limited observations for forecast verification, the seasonal forecast benchmarking and diagnosing task proves challenging. In this study, we develop a probabilistic deep neural network model, drawing on a wealth of existing climate simulations to enhance seasonal forecast capability and forecast diagnosis. By leveraging complex physical relationships encoded in climate simulations, our probabilistic forecast model demonstrates favorable deterministic and probabilistic skill compared to state-of-the-art dynamical forecast systems in quasi-global seasonal forecast of precipitation and near-surface temperature. We apply this probabilistic forecast methodology to quantify the impacts of initialization errors and model formulation deficiencies in a dynamical seasonal forecasting system. We introduce the saliency analysis approach to efficiently identify the key predictors that influence seasonal variability. Furthermore, by explicitly modeling uncertainty using variational Bayes, we give a more definitive answer to how the El Nino/Southern Oscillation, the dominant mode of seasonal variability, modulates global seasonal predictability.


Dynamical seasonal forecast and its three barriers
Global atmosphere-ocean-land coupled general circulation models (GCMs) that simulate the dynamics and interactions of various climate subsystems are the workhorses for predictions on weather to climate scales. Over the past seven decades, advances in computation, observation, modeling, and data assimilation have gradually extended the range of useful deterministic forecasts toward the predictability limit [1]. This limit of range, set by the chaotic nature of the atmosphere, is estimated to be around two weeks [2]. Moving beyond this range, subseasonal to seasonal forecasts seek prediction sources from low-frequency climate signals, which give rise to predictability by constraining the climate variability through the forecasting period.
By better characterizing the key prediction sources and their impacts, GCM-based seasonal forecasts have demonstrated markedly improved skill, offering crucial benefits to a wide range of societal sectors [3]. Despite the achievements, we highlight the following three barriers that hinder the improvement of dynamical seasonal forecast.
The first barrier is model initialization. Accurate initialization is a prerequisite for dynamical forecasts. However, existing initialization strategies often fail to constrain a model's state to faithfully match observations without introducing state shocks [4,5]. Furthermore, when compared to the atmosphere, the assimilation methods and observation networks for the ocean, cryosphere, and land are typically underdeveloped [6], inhibiting the full utilization of prediction sources from these climate subsystems.

arXiv:2010.14610v1 [physics.geo-ph] 27 Oct 2020
The second barrier is model formulation deficiencies. GCMs are simplified representations of the climate system. Their backbones are dynamical cores that compute the atmosphere and ocean dynamics in discrete geogrids. Besides the resolved dynamics, a suite of empirical parameterization schemes is used to account for the effects of unresolved processes. The considerable structural and parametric errors in parameterization can severely impair a model's forecast through multi-scale interactions.
The third barrier concerns the insufficient sampling of forecast spread. In dynamical forecasts, uncertainties in initialization and model formulation can quickly grow and yield distinct, but possible results. To estimate this forecast spread, we usually create forecast ensembles by sampling plausible initial states and plausible model formulations. A key question is: how many ensemble members are required to accurately infer the forecast spread? There has been a growing recognition that initial states accommodating different climate signals can yield distinct forecast spreads, resulting in variations of predictability [7,8] and requirement on ensemble size. Currently, we can not always afford large ensemble forecasts to accurately estimate forecast uncertainties. Also, the dependency of predictability on the initial climate state is not well-considered in deploying forecast ensembles or guiding forecast-related decision makings.
The difficulties in surmounting the outlined barriers have considerably slowed the pace toward more skillful seasonal forecasts. To accelerate forecast system development, we should set up reasonable forecast benchmarks [9], thereafter apply these benchmarks to clarify the limitations imposed by each aspect of the forecasting barriers [10]. Unfortunately, given the huge cost to initialize and run large-ensemble forecasts, and short history of climate observations for forecast verification, the seasonal forecast benchmarking and diagnosing task remains challenging.

Learning from climate simulations for seasonal forecast
Two recent studies suggest that we can turn to the data-rich climate simulation "model world" to seek implications for seasonal forecasts [11,12,13]. For instance, Ding et al. [12] drew analog from existing climate simulations that resemble observed ocean status as forecast ensembles, and found those analog forecasts match dynamical forecast systems in seasonal prediction for sea surface temperature and precipitation. This method could be interpreted as a k-nearest neighbour approach [14] from a machine learning viewpoint. In the second study, Ham et al. [13] trained a deep neural network to predict the El Niño/Southern Oscillation (ENSO) in climate simulations, and found that, with moderate modifications and using relatively few observational records, the neural network model applies well in real-world forecast of ENSO, outperforming dynamical forecasting for lead time of up to one and a half years.
The rationale behind these approaches is that, seasonal forecasts are conducted with identical or similar GCMs used in climate simulations. One can build analogical models using machine learning techniques to simulate the seasonal variability signal in climate simulations, and thereafter apply this GCM-revealed relationship for practical forecast.
While these attempts open the opportunity of employing large-scale statistical models (i.e., deep neural networks) to learn from climate simulation big data for seasonal forecasts, we notice the following three gaps along this novel line of research. First, existing research focuses on limited regions or phenomena, and have not fully exploited the wealth of climate simulation information to improve seasonal forecast for a broad range of regions. Second, forecast uncertainty, which is indispensable in seasonal forecast and crucial in supporting forecast-related decision making, has not been rigorously represented. Finally, any solid progress toward improving seasonal forecast should contribute to tackling the three aforementioned forecast barriers noted in Sect. 1a. However, it is unclear how the analog seasonal forecast practices inform the identification or conquering of these forecast barriers.
To address these deficiencies, we develop a probabilistic deep learning model that leverages the rich information from climate model simulations to enhance seasonal forecast capability and foster forecast diagnosis. Specifically, we summarize how key prediction sources regulate seasonal variability in climate simulations by learning a generative model that samples possible seasonal-ahead climate status conditioned on the current state of the predictors. By benchmarking dynamical seasonal forecasts with the proposed model, we demonstrate that current dynamical seasonal forecast systems have not fully exploited the forecast capability of existing GCMs. The impact of initialization error, GCM formulation deficiency, and internal climate variability are revealed. The results help answer the following questions: 1. Can we extend the data-driven seasonal forecasting paradigm to quasi-global scale?
2. How to represent uncertainty in a data-driven seasonal forecast?
3. How does this methodology inform potential improvement of dynamical seasonal forecasts?
We organize the rest of the paper as follows. Section 2 introduces the probabilistic modeling framework. Section 3 describes the data and experimental design. Section 4 compares our model with dynamical forecast systems in retrospective forecast experiments. We discuss the implications for diagnosing and understanding dynamical seasonal forecast in Section 5. Conclusions are drawn in Section 6.

Probabilistic modeling framework
The model we develop here is a customized Conditional Variational Auto-Encoder (CVAE) model ( [15,16,17,18], [19]). It is a conditional generative model that marries probabilistic graphical models with deep learning. Some of the key ideas here, such as variational inference and deep neural networks, may not be familiar to the present climate community, and relevant progress from the machine learning community is often built upon a series of recent work. Here, we provide a principle-based, step-by-step derivation of the probabilistic modeling framework. The specific model parametric form, architecture, and training strategy are thereafter explained.
Given climate simulations generated by GCM ensembles, we sample snapshots of predictor X and predictand Y pairs from these simulations. For instance, X could be the spatial distribution of upper ocean heat profile that has long memory and extensive impact on seasonal variability [13,20], and Y could be the corresponding spatial distribution of seasonal-ahead mean state of precipitation or near-surface temperature some months later. Based on the extracted (X, Y ) pairs, we build a probabilistic forecast model by approximating the conditional probability distribution of P (Y |X). This model has two sources of stochasticity. First, we focus on the impact of the considered predictor X, and discard explicit representation of less-predictable processes, such as the high frequency atmospheric component. Second, we inject model formulation noise by sampling climate simulation data generated by different GCMs.
We apply a parameterized function p θ (Y |X) to approximate P (Y |X). The learning process consist of optimizing θ such that p θ (Y |X) is close to P (Y |X) as revealed by our reference climate simulation data. While a suitable objective for optimizing θ is to maximize its conditional log likelihood log p θ (Y |X), a direct optimization is difficult given the high dimensionality and potentially arbitrary dimension dependencies of X and Y .
To alleviate this difficulty, we introduce a latent variable Z to govern the distribution of Y conditioned on X. This enables decomposing p θ (Y |X) into simpler building-blocks of p ψ (Y |X, Z) and p(Z|X): Here p ψ (Y |X, Z) is the conditional likelihood of Y given X and Z, with its parameter denoted by ψ; p(Z|X) is the conditional prior distribution of Z. We assume p(Z|X) to be standard normal distributed, and drop its dependency on X. This assumption is acceptable with flexible enough p ψ (Y |X, Z), which allows mapping the stand normal distributed Z to any preferred latent distribution form [21], meanwhile capturing the interaction between X and Z for inferring Y . Once p ψ (Y |X, Z) is estimated, we can generate samples of Y by sequentially sampling p(Z|X) and p ψ (Y |X, Z).
To apply maximum conditional likelihood estimation based on Equation 1, we need to integrate over Z, which is computationally intractable. To tackle this problem, we resort to stochastic gradient variational inference [15,17]. Specifically, we narrow the integration space of Z in Equation 1 by only considering Z values that are likely to generate Y given X. Such information is given by the unknown posterior distribution of P (Z|X, Y ). In variational inference, we bring in an inference probability function q φ (Z|X, Y ) with parameter φ to approximate P (Z|X, Y ). The error of this approximation is measured by the residual term as shown in the following equation (see appendix for derivation): Residual term (2) where E Z∼q φ [·] denotes expectation over q φ (Z|X, Y ), KL(· ·) is the non-negative Kullback-Leibler divergence of two distributions. Equation 2 stands at the core of variational inference theory. It decomposes the conditional likelihood into the reconstruction, regularization, and residual term, with the former two terms together coined as the ELBO. To train our model, we maximize the ELBO to implicitly maximize log p θ (Y |X) and minimize the residual term. The reason is two-fold. First, the ELBO stands for a lower bound estimate of log p θ (Y |X), with its tightness quantified by the residual term. By maximizing the ELBO, we will approximately maximize log p θ (Y |X) while strengthening the tightness of this lower bound estimate. The second reason is for computational convenience: by selecting proper distribution form and adopting the reparameterization trick [15,16], ELBO could be differentiable with respect to its parameters θ = {φ, ψ}, which allows scaling up the model to massive data via stochastic gradient optimization [22,23].
We illustrate the latter aspect by specifying the distribution form and the reparameterization trick. Following [15,16,17], we assume q φ (Z|X, Y ) to follow a multivariate normal distribution, with its mean vector and diagonal covariance matrix parameterized by two neural network functions, µ NN φ (X, Y ) and Σ NN φ (X, Y ). The neural network parametric forms are specified in the following section. With this assumption, the regularization term has the following analytical form, which can be optimized using stochastic gradient optimization: Here 1 is identity matrix, tr(·) and | · | are the trace and determinant operator, k is the dimensionality of the multivariate normal distributions.
To carry out gradient-based optimization for the reconstruction term, the reparameterization trick is applied to substitute the latent random variable Using this change-of-variable technique, we obtain a Monte Carlo estimate of the reconstruction term: Here denotes unbiased estimation, l is mini-batch size. Without loss of generality, we assume multivariate normal distribution of p ψ (Y |X, Z) = N µ NN ψ (X, Z), Σ ψ , where µ NN ψ is a neural network function that takes into input of X and Z to regress Y , Σ ψ is the covariance matrix. With this, we have: The right-hand side of Equation 5 is differentiable with respect to {φ, ψ}, allowing stochastic gradient minimization. It is worth noticing that Σ ψ balances the reconstruction term and the regularization term in composing the ELBO. Its impact is reflected in a hyperparameter β (see Table 1 in the latter context, [24]) in the final model.
To wrap up, Equations 2-5 offer an efficient variational inference framework by applying stochastic gradient ascent toward maximizing a lower bound of the conditional log-likelihood. The final loss function L of our learning model takes the following form: The first term encourages accurate reconstruction of the predictand using (1) the predictor information and (2) the latent variable that optimally explains the predictand. The second term forces the posterior of the latent to fully exploit the prior distribution space, encouraging each sample from the prior to be meaningful for generating plausible predictands. This gives the premise that, without knowing Y in prior, we can have good estimate of P (Y |X) by sequentially sampling the conditional prior of p(Z|X) and the likelihood of p ψ (Y |X, Z). The sampling results are thereafter used for probabilistic forecast.

Neural network parametric form
We specify the neural network parametric form of µ NN φ (X, Y ), Σ NN φ (X, Y ), and µ NN ψ (X, Z) here. Neural networks are powerful function approximators that apply multiple processing layers to learn hierarchical feature representations of data [25]. The computation units in a neural network can be flexibly organized to accommodate specific data structures and characteristics. Here we adopt convolutional neural networks (CNNs, [26]) as building blocks for our model. CNNs apply a hierarchical set of parameterized kernels to cross-correlate with the input feature maps to explicitly exploit the translational equivariance and compositionality properties of the data. Each cross-correlation is realized using a convolution operator: Here F L (x, y) and F L+1 (x, y) are the L-th and (L + 1)-th network layer feature vector at spatial coordinate (x, y), w L and b L are the weight and bias parameter of the L-th layer convolution kernel, whose receptive field size is m × n. Equation 7 enforces translational equivariance, which tells that, by moving every point in the input feature map by a given distance along a given direction (translation), the output feature map should in response move by a same distance along a same direction. This property is desirable as we prefer to maintain and take use of the spatial structure of the data through hierarchical feature learning, and efficiently recognize features shared across the space. Compositionality refers to the property that, complicated, large scale patterns are built from simpler, small scale ones. In CNNs, we learn compositional features by composing several convolutional and optionally down-sampling layers, obtaining a deep neural network that extracts local features in lower layers, while synthesizing local features for large-scale features in higher layers. The two peculiar characteristics of CNN can find analogues in numerical geofluid dynamics solvers [27], where we approximate the local non-linear behavior of fluid dynamics in a convolution manner, and apply numerical integration to obtain compositional representations of the fluid's spatiotemporal patterns. Recent years have witnessed increasing popularity of CNN applications in geophysical modeling [28,29,30,31,13].
Specific to our model, we apply CNNs to learn feature representations of X and Meanwhile, to reduce model complexity, we re-use the CNN extracted feature of X in µ NN ψ (X, Z), where the CNN extracted feature of X and the latent Z are concatenated and fed into a transposed convolution network to generate high-resolution predictand from low-resolution representations. The transpose convolution is implemented by first padding an input grid's neighbours with 0, followed by a regular convolution. This model architecture is illustrated in Figure 1.
It is worth noticing that, as we apply conventional CNNs on plane projections of spherical climate signals, we can not strictly guarantee translational equivariance, due to (1) distortion at polar regions, and (2) longitudinal discontinuity at the edge of the projection. Regarding the distortion problem, we do not consider the predictor/predictand information from polar region for the current work. This is acceptable, since most of the sources of seasonal predictability are believed to reside in the tropical oceans [32]. We alleviate the impact of discontinuity by rearranging the projection of the predictor to start from 20 • E, which is mostly masked by non-informative 0 values from the land areas, see predictor map in Figure 1.

Model training and forecast
We apply stochastic gradient descent [22] for model training. In preliminary experiments, we found that a joint training of µ NN φ , Σ NN φ , and µ NN ψ toward minimizing L occasionally yield poor prediction performance for the validation set. We attribute this to the fact that the model may fail to fully exploit the predictor information from X, but rely too much on the latent information Z to estimate Y . To have Z account only for the uncertainty information of Y conditioned on X, we develop a unique three-step training strategy. First, we set Z = 0, and carry out a preliminary training of µ NN ψ . This is equivalent to regressing Y with X. Second, we fix the network parameters in feature extraction of X, and update the rest of parameters of the model toward minimizing L. Finally, we fine-tune all the model parameters.
Note that the function of the encoder in our model is to find optimal transformation that maps the conditional prior of p(Z|X) to the unknown distribution of P (Y |X). This is achieved by (1) seeking corresponding sample space in the conditional prior of p(Z|X) that optimally explains the considered sample, and (2) fully exploiting the conditional prior to have each of the sample being meaningful in generating realistic Y samples. To make probabilistic predictions without knowing Y in prior, we first sample latent variable Z from conditional prior of p(Z|X), thereafter generate prediction Y base on p ψ (Y |X, Z). A common strategy is to generate multiple independent samples to obtain a probability distribution as prediction. Uncertainty is inferred from the forecast spread.

Experimental design and data
As a proof-of-concept, we consider a case study of seasonal forecasting for boreal winter (October-March) mean precipitation and near-surface (2m) air temperature. Following [13] and [20], we use the upper ocean potential temperature profile from the previous July as predictor. Both predictor and predictand are regridded to 2 • × 2 • , with a quasi-global coverage of 60 • S − 60 • N. Other predictors and predictands can be easily incorporated into the proposed framework.
We train our model using climate simulation data from GCM ensembles. We drive the trained model using ocean reanalysis for practical forecasts for the period 1982 to 2018. The forecasts are verified against observational precipitation or temperature records. We compare the performance of our model with calibrated dynamical forecasts using commonly-applied deterministic and probabilistic skill metrics. The applied datasets, dynamical forecast baselines, data processing methods, model architecture, and skill evaluation approaches are described below.
, takes into input of the latent variable Z and predictor X to generate samples of Y . L X and L Y denote the learned feature representations of the predictor X and predictand Y using convolutional neural networks. L X is shared by q φ (Z|X, Y ) and p ψ (Y |X, Z), as denoted by the black arrows. The blue/orange/red lines denote convolution/transposed convolution/dense layer operators in neural network modeling. The purple lines denote the reparameterization trick, showing how to generate posterior samples from q φ (Z|X, Y ) based on k-dimensional multivariate normal distributed samples and the inferred mean µ NN φ and covariance matrix Σ NN φ . Note that the encoder is only applied during model training.

Climate simulation
The climate simulation data for model training are obtained from the Phase 5 [33] and Phase 6 [34] of the Coupled Model Intercomparison Project (CMIP5/CMIP6). To achieve sample diversity and support the training of large-scale statistical models, climate simulations from various forcing scenarios are considered, including historical climate change simulations (in response to natural and human forcings), pre-industrial control simulations (representing natural climate variability that occurs in absence of external forcing agents), 1 percent/yr CO2 increase simulations, abrupt 4xCO2 simulations, and projections of future climate change guided by various representative greenhouse gas concentration pathways (for the CMIP5 projections) or shared socioeconomical pathways (for the CMIP6 projections).

Observational references
We use the ECMWF ocean reanalysis system 5 (ORAS5, [35]) July ocean potential temperature data as input to force our trained model for practical forecasts. To verify models' prediction skills, we compare models' forecasts with precipitation observational records from the Global Precipitation Climatology Project (GPCP) [36], and 2m air temperature reanalysis records from the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalyses (ERA5, [37]). The boreal winter mean observation are obtained by averaging the monthly data from October to March for each year from 1982 to 2018.

Dynamical forecasts
We compare our model's forecast performance with the North American Multi-model Ensemble (NMME) Phase-II retrospective forecast [38]. We consider four NMME dynamical seasonal forecast systems, namely the Community Climate System Model Version 4 (CCSM4, [39]) developed by National Center for Atmospheric Research, the third and fourth generation of Canadian Coupled Global Climate Model (CanCM3/CanCM4, [40]) developed by the Canadian Centre for Climate Modelling and Analysis, and the Forecast-oriented Low Ocean Resolution Model (FLOR-B01, [41]) developed by the Geophysical Fluid Dynamics Laboratory. Each model has 10 or more ensemble members. For consistency, we select only 10 ensemble members from each model. All forecasts have been bias corrected using cross-validation [42]. We consider the forecasts that are initialized in July from 1982 to 2018, and use their regridded October to March mean forecast as our baseline.

Data processing
Both predictors and predictands are spatial maps, regridded to a common 2 • × 2 • resolution using the nearest neighbor method. The predictor variable of ocean potential temperature at all depth layers is normalized using a common minmax normalization by subtracting the min and dividing by (max-min). The min and max scaling factors are set to 272.15 K and 310.15 K. Data for the land areas are masked by 0. The predictor is projected to start from 20 • E, which is mostly masked by non-informative 0 values from the land areas. The predictand variables are normalized at each grid point by subtracting the spatial mean, and dividing by the standard deviation of the observational records from 1982 to 2018 at grid scale.
We divide the climate simulation data into the training, validation, and test sets. The test sets are composed of the CMIP5 and CMIP6 historical climate simulations performed with each GCM from 1982 to 2007/2013. These data are selected to maximally match the practical forecast verification period. 90% and 10% of the rest data are shuffled and allocated to the training and validation sets.

Model architecture, hyperparameters and training
The structure, hyperparameters, and training details of the applied model is illustrated in Table 1. We use the Adaptive Moment Estimation (Adam, [43]) stochastic gradient descent for training. For each of the three steps in training, we set the maximum number of training epochs to 200, and employ an early stopping strategy: training is stopped when the relative decrease of the validation loss is less than 0.1% for 5 consecutive training epochs.

Skill metrics
We use the Pearson correlation coefficient score (r) and normalized root mean square error score (NRMSE) to measure the deterministic forecast skill of the ensemble mean predictions. In addition, we use the area under the ROC (Receiver Operating Characteristics) curve score (AUC, [45]) and the continuous ranked probability skill score (CRPSS, [46]) to assess the probabilistic skill of the full ensemble of forecasts. A summary of the skill metrics, their computation methods, and ranges is given in Table 2. † † L X and L Y are the convolutional neural network feature representations of predictor X and predictand Y . ‡ The selected hyperparameters are labeled with underline. ‡ ‡ D E and D D denote the depth of the convolution blocks. We sequentially double/halve the convolution channel in the encoder/decoder. § Each convolution block is composed of a convolutional layer, rectified linear activation (ReLU(x) = max(0, x)), batch-normalization [44], and 2 × 2 maximum pooling. § § β balances the reconstruction and regularization terms in the evidence lower bound equation.
Regarding the significance tests for the skill differences between models, we notice that the seasonal predictability exhibits strong annual variability, inhibiting the establishment of asymptotic distribution forms of models' skill score statistics. Thereafter, we apply bootstrap for significance test of models' forecast skill differences. We first generate 100 bootstrap samples of forecasting cases using random sampling with replacement. Thereafter, we calculate the ratio of forecast skill scores between two models for each bootstrap sample. The confidence interval is thereafter determined using a percentile method. All tests are single-side at confidence level of 95%.

Results
We first compare the performance of the NMME dynamical forecast systems and the CVAE model for predicting the observed October-March mean precipitation (Fig. 2) and 2m air temperature (Fig. 3). All forecasts are made from the previous July, with forecast lead time of 3-8 months. Evaluations are carried out at 2 • × 2 • grid scale for the period of 1982 to 2018. For the NMME forecasts, we consider single model 10-member based forecast, the intra-model mean skill of the four NMME models are shown in Fig. 2a,b and 3a, b. We also evaluate the four-model ensemble (4 × 10-member) forecast, results are shown in Fig. 2c, d and 3c, d. For the CVAE, we apply a 1000 ensemble member forecast created by sampling the latent variable Z. The results are shown in Fig. 2e, f and 3e, f. In these evaluations, we show the r skill score of the ensemble mean forecasts (left column of Fig. 2 and 3), as well as the AUC score that summarizes ensemble forecasts' hit rates and false alarm rates for a range of probability thresholds (right column of Fig. 2 and 3). Stippling denotes grids where the NMME and CVAE forecasts show significant skill difference at 95% confidence level. The spatial average skill score and the number of stippling gridboxes are used to roughly Step 1: Construct the ROC ‡ curve by (0.5,1) Area under the curve drawing hit rate against false alarm rate (AUC) for predictions exceeding the ten deciles.
Step 2: Compute area under the curve. Continuous ranked probability skill score (0,1) (CRPSS) † Min/max value indicates no skill/perfect skill. † † n is sample size, y i simu , y i obser are individual simulation/observation samples indexed by i. Overbar denotes climatology mean. ‡ Receiver Operating Characteristics. See appendix for detailed explanation. ‡ ‡ P simu /P obser is the cumulative probability distribution function of the simulation/observation. P obser (y) = 0, if y < y obser 1, otherwise . P ref is the cumulative probability distribution of climatology records. See [7] for computation details.
summarize models' forecast capability. Evaluation using NRMSE and CRPSS shows similar patterns, and is given in the Appendix.
Compared to the NMME dynamical forecasts, the CVAE model exhibits similar skill distribution pattern, but achieves an overall higher level of deterministic and probabilistic skills. Specific for precipitation forecast (Fig. 2), the CVAE forecasts show significantly higher r and AUC skill score in Eurasia, Eastern and Southern Africa, tropical Pacific, Southern U.S., and Amazon (Fig. 2e, f), compared to either the intra-model mean skill score of the NMME forecasts (Fig. 2a, b), or the skill score of the NMME multi-model ensemble forecast (Fig. 2c, d). While the multi-model ensemble can considerably enhance the r skill score of dynamical forecasts (Fig. 2a, c), reducing its gap with CVAE forecasts (Fig. 2e), we notice that the AUC scores of the multi-model ensemble is worse than individual models for most regions, which may due to the climatology mismatch of the considered dynamical forecast models.
For 2m air temperature forecast (Fig. 3), both NMME and CVAE forecasts achieve higher skill scores compared to precipitation forecasts. The CVAE model demonstrates particularly strong skill in predicting the temperature variation for the eastern Pacific ENSO region, with r > 0.9 (Fig. 3e) and AUC> 0.85 (Fig. 3f). Also, the CVAE model shows significantly higher prediction skill for regions where the NMME skill is limited, such as western Eurasia and central North America.
It is worth noticing that the skills of the CVAE model shown here are derived purely from GCM simulations, no observational data has been used for training the CVAE model, ensuring that we run a small risk of "overfitting", in the sense that the CVAE predictions reflect quirks of the climate simulation data, rather than regularities generalizable to the real-world climate system. Based on these results, we have shown that the CVAE model offers a competitive benchmark for both deterministic and probabilistic seasonal forecasting.

Forecast diagnosis
How can we outperform GCM-based dynamical seasonal forecast using GCM simulation information only? What are the implications for diagnosing and improving dynamical seasonal forecast? We answer these questions by contrasting a same-GCM supported dynamical forecast and CVAE forecast, applied to either the real-world climate system or the GCM simulation "model world". Corroborated by multi-GCM supported forecasting results in Fig. 2 and Fig. 3, we reveal the forecasting barriers posed by GCM initialization and formulation deficiencies, as well as this GCM's predictability limit.
We choose the Canadian Earth System Model (CanESM) and its corresponding seasonal forecast system CanCM4 as an example. We consider the following two models: (1) the CanCM4 dynamical forecast system, and (2) a new  . (a, b), the intra-model mean r and AUC skill score of four NMME models. (c, d), the r and AUC skill score of the NMME four-model ensemble.
(e, f), the r and AUC skill score of the CVAE model. Stippling in (a-d) shows locations where the NMME forecasts show significant (95% confidence level, same for the rest tests) higher skill than CVAE forecasts. Stippling in (e) shows locations where the CVAE forecasts show significant higher r score than NMME ensemble forecasts.
Stippling in (f) shows locations where the CVAE forecasts show significant higher AUC score compared to the mean AUC score of four NMME models. The spatial average skill score and the number of stippling are denoted by µ and n in each sub-figure. CVAE model trained with the 12,253 years of the CanESM simulation data (CVAE CanESM ). We apply CVAE CanESM for forecast in (1) the observed real-world climate system, and (2) the "model world" realized by the CanESMv5 historical simulation. These forecasts are denoted as CVAE Obser CanESM and CVAE CanESMhis CanESM , respectively. Results measured using r skill score are shown in Figure 4. Comparing the CanCM4 and CVAE Obser CanESM highlights potential initialization errors, while contrasting CVAE Obser CanESM , CVAE CanESMhis CanESM and CVAE Obser Ensemble (our CVAE model based on multiple CMIP5 and CMIP6 simulations) highlights the model formulation deficiency and predictability limits.

Initialization error
Since both the CanCM4 and the CVAE CanESM are based on the same GCM, their performance difference in practical forecast should largely depend upon how they make use of this GCM. In CanCM4 forecasts, we run simulations starting from an ensemble of initial states that agree with observational constraints. In CVAE Obser CanESM forecasts, we make use of the seasonal variability signal learned from continuous runs of CanESM simulations. Although the probabilistic relationship learned by the CVAE CanESM may be oversimplified or misrepresented, the benefit of this simplification lies in that, it allows bypassing iterative state updating in forecast, thereafter greatly alleviate the drift problem in initializing dynamical forecasts.
The discussion above tells that, the skill difference between the CanCM4 and the CVAE Obser CanESM reflects the following tug-of-war: Does the initialization error in dynamical forecasts surpass the error of the CVAE model resulting from oversimplifying or misrepresenting the dominant factor of seasonal variability? For precipitation forecasts (Fig. 4a,  c), while CanCM4 demonstrates significant better r skill score for certain regions, such as tropical and western mid- Figure 4: Correlation coefficient (r) skill score of CanCM4 and CanESM-supported CVAE forecasts. (a, b), r skill score for CanCM4 precipitation/temperature forecast. (c, d), r skill score for CanESM-supported CVAE precipitation/temperature forecast. (e, f), r skill score for CanESM-supported CVAE forecast applied to CanESM historical simulations. Stippling in (a, b) denotes locations where the CanCM4 forecasts show significant (95% confidence level, same for the rest tests) higher skill than corresponding CVAE forecasts. Stippling in (c, d) denotes locations where the CVAE forecasts show significant higher skill than CanCM4 forecasts. Stippling in (e, f) denotes locations where the CVAE forecasts applied to the CanESM historical simulation show significant higher skill than the CVAE forecasts applied to real-world climate system. The spatial average skill score and the number of stippling are labeled.
latitude Pacific, there are many other regions where CVAE Obser CanESM obtains significant better r skill score, such as northern Eurasia, central U.S., and Southern Hemisphere mid latitude. For 2m air temperature forecasts (Fig. 4b, d), compared to CanCM4, CVAE Obser CanESM achieves significant better r skill score for a broader range of regions, including both continental and ocean parts. These results indicate that CanCM4 has not fully exploited the seasonal prediction capacity of its GCM. We believe this skill lagging of the CanCM4 compared to CVAE Obser CanESM offers a lower-bound estimate of unrealized prediction capacity of the GCM due to initialization deficiencies.

Model formulation deficiency and predictability limit
It has become a conventional practice to reveal the defect of a single GCM by comparing its forecast with forecast from multi-GCM ensembles. In dynamical forecast, this comparison reflects a compound impact of initialization and GCM formulation deficiencies. The previous section has demonstrated that the initialization deficiencies can be largely alleviated using the CVAE methodology. Regarding the GCM formulation deficiencies, it is beneficial to further clarify whether they originate from a biased representation of the seasonal variability signal, or an underestimation of this signal. We shed light on this problem by contrasting the following three CVAE forecasts: (1) CVAE Obser Ensemble , the multi-GCM ensemble based CVAE forecast applied to the observed real world (Fig. 2e and 3e), (2) CVAE Obser CanESM , the single GCM (CanESM) based CVAE forecast applied to the observed real world (Fig. 4c,

d), and (3) CVAE CanESMhis
CanESM , the single GCM based CVAE forecast applied to this GCM's historical simulation "model world" (Fig. 4e, f). The results reveal the achieved and potentially achievable seasonal forecast skill of the considered GCM.

While the CVAE Obser
Ensemble and CVAE Obser CanESM forecasts show considerably similar spatial distribution of skill ( Fig. 2e and Fig. 4c, Fig. 3e and Fig. 4d), we focus on the hot spots where these two forecasts show significant skill differences. Figure 5 shows the bootstrap distribution of spatial average r skill scores for four of these hot spots, namely Eastern Africa, Southeast China, Maritime Continent, and Eastern U.S. We distinct three configuration cases based on the ranking of the spatial average r skill score achieved by these CVAE forecasts.

CVAE Obser
Ensemble > CVAE CanESMhis CanESM > CVAE Obser CanESM . This ranking of r skill score (red rectangle delineated box charts in Fig. 5) suggests that, compared to the real-world climate system, the CanESM model alone fails to develop a strong enough constraint posed by the prediction sources in its simulations. This happens for precipitation forecast in Eastern Africa and Southeast China. For both cases, this overly weak constraint of CanESM is meanwhile biased when applied for practical forecast, given the negative r skill score of the CVAE Obser CanESM forecast.

> CVAE Obser
Ensemble > CVAE Obser CanESM . This ranking of r skill score (blue rectangle delineated box charts in Fig. 5) represents the most common case. It suggests that, compared to the real-world climate system, the CanESM model develops a stronger constraint posed by the prediction sources in its simulations. However, this strong seasonal variability signal does not apply well to the real-world climate system, compared to the seasonal variability signal revealed by the multi-GCM ensemble.

CVAE CanESMhis
CanESM > CVAE Obser CanESM > CVAE Obser Ensemble . This ranking of r skill score (green rectangle delineated box charts in Fig. 5) tells that, compared to the multi-GCM ensemble, the considered CanESM is advantageous in representing the seasonal variability signal, and achieves a better forecast skill for the selected region.
The discussion above suggests that, for the considered CanESM example, there are regions where the GCM develops overly weak and biased constraints on precipitation variability (red rectangle delineated box charts in Fig. 5a). We believe these are hot spots worth close and urgent attention in model development. Besides these regions, the prediction sources from upper ocean thermal profile usually pose a stronger constraint on seasonal variability in the GCM simulation than in the real-world climate system, especially for temperature ( Fig. 4f and Fig. 5b). This strong signal reflects the potentially achievable forecast skill of the considered GCM. While this predictability signal may vary with both externally forcings and internal climate variability, we examine how this signal depends upon the dominant seasonal variability mode of El Niño/Southern Oscillation (ENSO) in the latter section, and leave the discussion for the impact of externally forcings in future works. Overall, the comparison results here help to reveal the GCM formulation deficiencies, and pinpoint the directions toward forecast improvement.

Identifying sources of predictability
Identifying the sources of predictability helps specifying the key climate processes that offer predictability. In dynamical simulation, this task often requires costly computation as we examine GCM's responses to changes from many potentially crucial factors. Here, we alleviate this computation burden by analyzing the CVAE model instead. The analysis in turn tells if the CVAE model applies physically-meaningful factors in its forecast.
The technique we apply here is called saliency map analysis, which shows how each portion of the predictor influences the forecast by visualizing the neural network gradient. In other words, it is a visualization tool highlighting the geographical aeras that have influenced the forecast the most. Consider an example case of precipitation forecast in the Maritime Continent, we calculate the change of CVAE ensemble-mean precipitation forecast for this region due to a one standard deviation increase of potential temperature during the previous July at each location and ocean depth. Results for all the forecasting years , the strong El Niño years (Niño3.4 index in July is larger than 1), and the strong La Niña years (Niño3.4 index in July is smaller than -1) are averaged and plotted for three representative ocean depth layers (5m, 45m, and 150m) in Fig. 6.
The zonal dipole patterns in Fig. 6 suggest that, the CVAE precipitation forecast in the Maritime Continent is negatively influenced by the July east-central Equatorial Pacific top-layer ocean potential temperature, and positively influenced by the July local top-layer ocean potential temperature. These patterns faithfully reflect the dominant role of the Walker circulation. A strong Walker circulation, arising from an enhanced ocean thermal gradient in the  Pacific, shifts its rising branch to the east, and introduces anomalous sinking over the Maritime Continent, suppressing precipitation formation.
The well-pronounced impact at 45m ocean depth (Fig. 6d-f) highlights the role of below-surface ocean thermal structure for leading and sustaining the sea surface and atmosphere circulation variations at seasonal lead time. Furthermore, the comparison between saliency responses in El Niño years (Fig. 6b, e, h) and La Niña years (Fig. 6c, f, i) illustrates the asymmetric impact of ENSO, namely precipitation during El Niño is more sensitive to ocean thermal variations in the eastern Niño 3 region (grey rectangle box), especially for the deeper ocean, while precipitation during La Niña is more sensitive to ocean thermal variations in the central Niño 4 region (black rectangle box) and less influenced by the deeper ocean. These findings, corroborated by existing observation and simulation studies [47], confirm that the CVAE model leverages physically-meaningful patterns in its predictions. This saliency map analysis enables efficient identification of key prediction sources.

Impact of ENSO on predictability
Seasonal forecasts starting from different climate initial states may demonstrate distinct levels of skill, which has been partially attributed to the impact of internal climate variability on predictability [48,7,8,49]. For instance, some climate conditions may impose persistent constraints on seasonal variability, offering opportunities of high predictability; while some other conditions may amplify the state estimation uncertainty, and introduce forecast busts [50]. Retrospective forecasts employing limited ensemble members often fail to give conclusive answer to this predictability-state dependency relationship. Here, we examine if the CVAE model can leverage the rich information from climate simulations to shed light on this problem. The results in turn verify if we make reasonable forecast uncertainty estimations using the CVAE methodology.
The analysis here focuses on the dominant seasonal variability mode of ENSO and its impact on precipitation predictability. To infer this impact, we apply the forecast intra-ensemble variance as an indicator of predictability. We correlate this intra-ensemble variance with the forecast starting time ENSO status (measured using the Niño3.4 index) for the period of 1982 to 2017. Results are shown in Fig. 7.
The intra-ensemble variance reflects the growth of state estimation uncertainty in model state space. In NMME forecast, we explicitly compute this growth of state estimation uncertainty, starting from an ensemble of initial state estimates that are loosely constrained by observations. In CVAE forecast, we adopt a top-down strategy, that is, we discard explicit modeling of the processes that are less predictable or less relevant in seasonal forecasts, and represent their impacts using variational inference. Given their distinct representation of uncertainty, it is quite surprising that, despite the magnitude of difference, the NMME and the CVAE model show very similar correlation patterns between Niño3.4 and their forecast intra-ensemble variances (Fig. 7a, b). The relatively weak correlation between the NMME intra-ensemble variance and the Niño3.4 index (Fig. 7a) may reflect the ensemble size limitation in dynamical forecast. Based on limited ensemble members, the intra-ensemble variance estimation can be seriously influenced by sampling variability noise. To support this argument, we draw analog from the CVAE forecast: as we decrease the ensemble size of CVAE forecast from 1,000 to 10, the resulting correlation map (Fig. 7c) shows similar level of correlation magnitude as in Fig. 7a.
The discussion above suggests the necessity for applying large ensembles to infer the predictability limit and its connections with climate variability signals. While this practice is prohibitively costly in dynamical forecast, the CVAE model offers useful analogical insights. Compared to the NMME result, the correlation map between the CVAE intra-ensemble variance and the Niño3.4 index gives a more definitive answer to how ENSO modulates predictability. In Fig. 7b, for the negative correlated regions, such as Southeast Asia, Australia, and South Africa, forecasts starting in El Niño Julys tend to give a narrower ensemble spread in predicting the mean boreal winter precipitation, which offer an opportunity of predictability; contrarily, predictions starting in La Niña Julys tend to give a wider forecast spread, highlighting the necessity of having larger ensembles to sample the less-predictable climate. The opposite analysis applies for the positive correlated regions: El Niño July forecasts tend to give a wider ensemble spread and less predictability, while La Niña July forecasts tend to give a narrower ensemble spread and more predictability.
Finally, we notice a strong resemblance between correlating Niño3.4 with intra-ensemble variance (Fig. 7a, b, c) and correlating Niño34 with mean precipitation amount (Fig. 7d). This suggests a global pattern for the considered seasonal forecast setting: The abnormally high precipitation associated with ENSO is accompanied by a larger forecast uncertainty, while abnormally low precipitation associated with ENSO is accompanied by a smaller forecast uncertainty.

Concluding Remarks
We have developed a deep generative model that draws on a wealth of existing climate simulations for probabilistic seasonal forecasts. The model offers a competitive quasi-global scale seasonal forecast benchmark, and a powerful tool for disentangling the dynamical forecast limitations imposed by initialization errors, model formulation errors, and internal climate variability. Three of the key contributions here are summarized as follows: Derivation of the evidence lower bound equation Equation 8 gives a step-by-step derivation for the evidence lower bound equation. Compared to a commonly-used strategy of starting from the viewpoint of minimizing the Kullback-Leibler divergence between posterior P (Z|X, Y ) and the variational distribution families, the derivation here focuses on decomposing the conditional marginal log likelihood log p θ (Y |X) into the ELBO and the residual term. As the result suggests, by maximizing the ELBO, we will approximately maximize log p θ (Y |X) while minimizing KL q φ (Z|X, Y ) P (Z|X, Y ) .

Regularization term
Evidence lower bound Residual term (8) Climate simulation data

ROC curve and Area under the ROC curve score
To construct the ROC curve, we consider probability forecasts for the categorical events that the predictand exceeds pre-defined thresholds. The hate rate and false alarm rate are estimated for these events. For instance, consider the event that the predictand is larger than 20% quantile of the observations or forecasts, the hit rate (HR) and false alarm rate (FAR) are calculated as follows: HR 20% = n (y simu > y 20% simu ) ∧ (y obser > y 20% obser ) n(y obser > y 20% obser ) FAR 20% = n (y simu > y 20% simu ) ∧ (y obser ≤ y 20% obser ) n(y obser ≤ y 20% obser ) Here HR 20% and FAR 20% are hit rate and false alarm rate for the event that the predictand Y exceeds the 20% quantile, n(·) denotes the number of events, y simu > y 20% simu represents the event that the prediction is larger than the 20% quantile of all the prediction records. y obser > y 20% obser /y obser ≤ y 20% obser represents the event that the observation is larger/no larger than the 20% quantile of the observation records.
In Fig. A1, we show how to construct the ROC curve and calculate the AUC score using example of winter precipitation forecast for the Southeastern U.S. (23 • N-27 • N, 81 • W-97 • W). For each of the considered model, we calculate the (HR, FAR) for the events that predictand exceeds the 20%, 40%, 60%, 80% quantile of the observations or predictions. We join these points to obtain the ROC curve. The area under the ROC curve is thereafter applied to summarize how reliable the model is at predicting these categorical events. Seasonal prediction skill evaluation using supplementary skill metrics Figure A2 and A3 compare the NMME and CVAE forecast performance using supplementary deterministic and probabilistic skill metrics of normalized root mean square error (NRMSE) and continuous ranked probability skill score (CRPSS). Figure 9: Forecast skill for October-March mean precipitation using models starting in July for 1982 to 2018 at 2 • × 2 • grid scale, measured using supplementary skill metrics of normalized root mean squre error score (NRMSE) and continuous ranked probabilistic skill score (CRPSS). (a, b), the intra-model mean NRMSE and CRPSS skill score of four NMME models. (c, d), the NRMSE and CRPSS skill score of the NMME four-model ensemble. (e, f), the NRMSE and CRPSS skill score of the CVAE model. Stippling in (a-d) shows locations where the NMME forecasts show significant (95% confidence level, same for the rest tests) higher skill than CVAE forecasts. Stippling in (e) shows locations where the CVAE forecasts show significant higher r score than NMME ensemble forecasts. Stippling in (f) shows locations where the CVAE forecasts show significant higher AUC score compared to the mean AUC score of four NMME models. The spatial average skill score and the number of stippling are denoted by µ and n in each sub-figure.