Rare Event Simulation of Extreme European Winter Rainfall in an Intermediate Complexity Climate Model

We test the application of a rare event simulation (RES) algorithm to accelerate the sampling of extreme winter rainfall over Europe in a climate model. The genealogical particle analysis algorithm, an ensemble method that interrupts the simulation at intermediate times to clone realizations in which an extreme event is developing, is applied to the intermediate complexity general circulation model PlaSim. We show that the algorithm strongly reduces the numerical effort required to estimate probabilities of extremes, demonstrating the potential of RES of seasonal precipitation extremes.

2 of 11 are required to sample extreme seasons. At the same time, higher model resolution benefits the representation of storm tracks, storm structure, the water cycle and extreme precipitation, and a global coupled modeling setup allows for the representation of teleconnections with the tropics or the Arctic (M. D. Priestley & Catto, 2022;Roberts et al., 2018;Schiemann et al., 2018;Vannière et al., 2019;Zappa et al., 2013). Limited computing resource may require compromises between these different demands.
Efficient simulation strategies offer a method to resolve such dilemmas. Rare event simulation (RES) algorithms, as applied in this article, provide ways to generate many more samples of an extreme event. This larger set of samples in turn allows for more accurate estimates of the probability of occurrence of the event, or to obtain the same accuracy with less numerical effort. Originally developed in statistical physics (Kahn & Harris, 1951), RES algorithms have found applications in, for example, queuing networks (Heidelberger, 1995), biochemistry (Kuwahara & Mura, 2008), and aircraft design (Morio & Balesdent, 2015). In Earth system modeling, rare event algorithms have previously been successfully applied to atmospheric models, to over-sample heat waves (Ragone & Bouchet, 2021;Ragone et al., 2017) and extremely intense tropical cyclones (Webber et al., 2019).
The aim of this study is to demonstrate the potential of RES for simulating seasonal precipitation extremes. While RES has already been successfully demonstrated on Earth system models before, different types of extremes develop differently over time. A successful RES setup for one quantity, such as temperature, may not necessarily carry over, with the same efficiency improvements, to other quantities, such as precipitation. As experiments with high resolution climate models are numerically demanding, there is a need for exploratory demonstrations of RES on intermediate complexity models. We therefore implement RES of winter rainfall in an intermediate complexity climate model, as a test bed for further studies on higher resolution models. We will estimate probabilities of threshold exceedance of winter precipitation far into the tail of the distribution and quantify how our RES approach has significantly reduced the numerical effort required, compared to traditional brute force approaches.
We will use an ensemble based genealogical algorithm similar to those used in Ragone and Bouchet (2021), Ragone et al. (2017), andWebber et al. (2019), which are, in turn, variants of the algorithm described in Del Moral and Garnier (2005). In this family of algorithms, an ensemble is run until an intermediate stopping time, at which time an objective function is applied to each realization. Realizations with a higher value of the objective function are deemed to be more promising and are cloned. Realizations with low values may be terminated, hence focusing numerical effort on the extreme event.
In Section 2, we will describe the atmospheric model used, introduce rare event simulations methods in general and the setup used in this study in particular. In Section 3, we will then present the outcome of the experiment, demonstrating that the method used can generate a much larger sample of extreme events, reducing strongly the total amount of simulation time needed to obtain the same accuracy, as quantified by the relative sampling error of the estimate.

The General Circulation Model
The model used in this study is the open-source Planet Simulator (PlaSim) model (Fraedrich et al., 2005), which uses the Portable University Model of the Atmosphere as its dynamical core. PlaSim gives a reasonable representation of atmospheric dynamics and of their interactions with the land surface and with the mixed layer of the ocean; it includes parameterizations of radiative transfers, diagnostic cloud cover, large-scale precipitation, convective precipitation, and dry convective adjustment.
The model is used in a T42 horizontal resolution with 10 vertical layers, fixed default model parameters and climatological SST and sea ice extent with a seasonal cycle. The time step is 20 min.
The model climate has been shown to represent basic features of the general circulation, including the North Atlantic storm track (Andres & Tarasov, 2019) (see also Figure 6).

Importance Sampling
The fundamental idea behind many RES methods is that of importance sampling (IS). Assume we wish to estimate a small probability = ℙ{ ∈ } for a random variable X. The set ℝ could be, for example, those values exceeding a high threshold A = (θ, +∞), or a small interval of width Δθ around a target value A = (θ − Δθ/2, θ + Δθ/2). The random variable X could be a physical observable, such as temperature, wind intensity or precipitation, possibly geographically, and temporally averaged.
A method to estimate γ is to generate N independent samples X i of X and use the brute force Monte Carlo estimator where is the indicator function on the set A ( ( ) = 1 for x ∈ A and 0 otherwise). The problem with this approach is that the relative error (RE), defined as REMC = √ (̂)∕ , the standard deviation of the estimator relative to γ A , is approximately 1∕ √ for small γ A (Rubinstein & Kroese, 2017). Therefore, for rare events with γ A ≪ 1 we would need a very large number of samples N to keep the RE at an acceptably low level.
A solution is to sample a different random variable ̃ whose distribution is linked to that of X in a known relation. Let's assume that both X and ̃ have a continuous distribution with probability density ρ(x) and ( ) respectively, with such that ( ) > 0 whenever ( ) ( ) > 0 . By generating independent samples ̃ from ̃ the IS estimator̃= can be implemented, where ( ) ∶= ( )∕̃( ) whenever ( ) > 0 and zero otherwise. It is important to note here that the densities ρ and do not need to be known. Only the proportionality factor L(x) needs to be known, which, as we will see below, is the case in the algorithm that we use.
The IS estimator is still unbiased and by choosing a density that makes the rare event A more common, the RE of this estimator, REIS = √ (̃)∕ , can be greatly reduced in comparison to that of the brute force Monte Carlo estimator (RE MC ). The exact reduction depends on the probability distributions involved and can rarely be computed analytically. At least theoretically, an IS estimator with zero RE can be constructed by taking = ( ∕ ) , although it is of no practical use, since it involves knowledge of γ A , the very quantity we want to estimate. This optimal change of measure does illustrate that the change of measure should make the rare event more common. It should also be noted that a poor choice of distribution (one that makes the rare event even rarer), can lead to an increase in RE. For more details on IS, see, for example (Bucklew, 2004;Rubino & Tuffin, 2009;Rubinstein & Kroese, 2017).

Rare Event Simulation Algorithm
When running a brute force simulation of a dynamical process, like the Earth's climate system, the model spends most of the time simulating a completely normal state of the climate, perhaps 99.9% of the time when considering a 1/1000 years event. The aim of RES algorithms is to perform the change of measure that is necessary to implement IS, spending more effort simulating relevant realizations. The random variable X is now a random process in time and the probability distributions ρ and should be thought of as probability distributions on possible paths of this process.
The idea of IS can be generalized to stochastic processes. In such cases, a change of measure can be carried out by altering the dynamical equations defining the process, for example, by changing the drift of a diffusion process. The probabilities of sample paths is then altered, with the Radon-Nikodym derivative given by the Girsanov theorem. This can again be used to implement an unbiased estimator, as in the univariate case above (Vanden-Eijnden & Weare, 2012). When considering a high-dimensional climate model, this method faces several challenges. First of all, it is often not evident how the system should be altered to generate more realizations of the rare event Instead, we can run an ensemble simulation with a number of parallel model runs, and at intermediate times interrupt the calculation to perform a selection procedure. This way, some paths will become more likely than others, performing a change of probability measure. Since we can keep track of how much more likely the event has become due to the algorithm, we can find the proportionality factor L, necessary to construct an unbiased estimator from the manipulated ensemble.
In the class of algorithms known as genealogical particle analysis (GPA) (Del Moral & Garnier, 2005), trajectories that are not developing a rare event have a higher likelihood of being terminated. In our application this would be ensemble members with little precipitation over Europe. The most promising trajectories are, on the contrary, cloned, with more clones assigned to realizations with more precipitation. The algorithm is presented in Appendix B.
The central idea behind GPA is the following: most of the time the different realizations follow a free evolution (step 2a in Appendix B). At intermediate selection times t k , the ensemble is manipulated to make the rare event of interest more common. At step 2b a weight is calculated, which determines how many clones will be created in step 2c. The number of clones can be zero, meaning this realization will be terminated. Note that no changes to the model code are necessary to implement GPA, which is an important characteristic when working with complex numerical codes.
The selection steps 2b and 2c are where the change of measure is performed. As in the standard IS estimator 2, the change of measure needs to be taken into account in the estimation step 4. The probability of each realization has been changed at each selection step by a factor ̄( ) , so the overall probability has been changed by the product This, therefore, needs to be divided out to make the estimation unbiased.
Importantly, an estimator of the variance of the GPA estimator, and hence the RE, has been derived (Del Moral & Garnier, 2005), allowing us to compare a posteriori the accuracy of the GPA estimator to brute force Monte Carlo (see Equation B2 in Appendix B).
For extremes of time-integrated observables ∫ 1 0 ( ( ))d , such as the time-integrated precipitation, an exponential weight function can be used, favoring realizations with a sustained high value of the integral of V. The control parameter C determines the strength of the selection procedure, allowing the user to focus on more or less rare events. Figure 1 illustrates the algorithm. As the algorithm doesn't require any modification of the climate model, only terminating and cloning of realizations, it can be applied to complex climate models.

The RES Configuration Used
The random variable we are interested in is the seasonally integrated total precipitation over a given area is the state of the model at time t, p(x(t), r) is the precipitation rate at geographical position r at time t, Ω ⊂ S 2 is the geographical region of interest and T 0 and T 1 are the initial and final times, respectively. Here Ω is taken to be the region between −20°E and 50°W and 30° and 70°N. The temporal integration is performed over December, January and February, that is, between T 0 = 1 December 00:00 and T 1 = 1 March 00:00.
The initial distribution is sampled M = 250 times on 1 December from a 500 years control simulation. Cloning is done every 5 days, that is, t k − t k−1 = 5 days. This cloning interval is chosen such that clones have sufficient time to diverge in between cloning, but don't revert to the mean. The error doubling time for PlaSim is estimated to be 2.7 days (Lyu et al., 2018). As control simulation we perform a stationary run (with seasonal cycle), using the default parameter settings of the PlaSim model.
For the weight function W k we use the same observable as the one for which we are estimating rare events, that is, The parameter C was chosen based by fitting a Gamma distribution to a histogram of the variable of interest , total winter precipitation over Europe, in the control run and picking the value of C that would result in a change in the mean of the RES ensemble of 3 standard deviations (see Appendix A). The assumption of a Gamma distribution is never used in the algorithm, it is only used to obtain a reasonable value of C to target a given range of extreme events. The resulting value is C = 371.30.

Effectiveness of the RES Algorithm
We first verify that the RES is working effectively. Figure 2a shows the cumulative rainfall anomalies generated in the RES versus the brute force control run. We see that, as expected, over time the realizations with the highest precipitation have been selected, tilting the distribution at the final time toward extreme events. Figure 2b shows histograms of winter precipitation generated from 490 years of the 500 years control run (discarding the initial 10 years) and from RES. It is evident that the RES has sampled the tail of the winter European rainfall distribution. At the final time what was a rare event in the control run has become common in the RES. The mean rainfall level in the RES ensemble lies at 0.220 m. Only 13 winters in the 490 years control run have winter rainfall exceeding this level.
Using the ensemble generated in the RES run we can estimate probabilities using the estimator B1 and compare them to those obtained from the brute force simulation using estimator 1. Figure 3 shows the estimated complementary cumulative distribution function ̄( ) = ℙ(̄) , giving the probability of exceeding a threshold θ. We see that for all estimated probabilities, the estimate for the RES is within the ±2σ interval of the brute force simulation, even though the RES simulation has a total simulation time of only 61.45 years, compared to the much longer control run of 490 years. Figure 4 shows estimates of the normalized RE for the brute force estimator and for the RES. The normalized RE is the RE for sample size N = 1, hence comparing relative errors for equal sample sizes. Figure 4 demonstrates that in the range that the RES is preferentially sampling (roughly 0.21-0.23 m, see Figure 2b) there is a strong reduction in the normalized RE. This reduction demonstrates that with the same numerical effort, we can get a much more accurate estimate using RES, or a similar accuracy with less numerical effort. This is the case, at least, in the range the RES is preferentially sampling. By varying the parameter C, this range can be tuned to a range of interest. As the RE of the brute force estimator is proportional to 1∕ √ , an increase of N to r 2 N would be required to reduce the RE by a factor r with brute force Monte Carlo. Said otherwise, if RES can reduce the normalized RE by a factor r, this translates into a reduction of computational effort by a factor r 2 . For the rarest events (in the range (0.2275, 0.2325)), RE is reduced by a factor 7.53, meaning the numerical effort could be reduced by a factor 56.84.
A complementary view to Figure 3 is presented in Figure 5. Here, interval probabilities ℙ( − Δ ∕2 <̄< + Δ ∕2) are estimated. We see that both estimates fall within each other's 2σ interval. As expected from 4, the width of the 2σ interval increases faster for the brute force estimates than for the RES estimates. Specifically, we can see from Figure 5 that, for example, the lower bound of the 2σ interval for the estimated probability for the interval (0.225, 0.230) becomes negative for the brute-force simulation, while the RES still gives a positive lower bound of 0.00559. By definition the probability estimates of the brute-force MC are positive (see Equation 1). The negative lower bound of the 2σ interval is a manifestation of the increasing RE caused by a high probability that, at this sample size, an event of this rarity will not be realized, resulting 6 of 11 in a zero probability estimate. It should be noted that this occurs when the inverse probability becomes of a similar magnitude to the sample size.

Characterization of Extreme Winter Rainfall Condition
An important advantage of RES methods, compared to extrapolations based on extreme value theory, is the availability of the full dynamical fields resulting from the numerical integration of the RES ensemble. Here we illustrate this property by comparing the physical fields of the RES ensemble to those of the control run. Figure 6 shows the mean circulation and precipitation over the North Atlantic and Europe as represented by the PlaSim CONTROL and RES experiments. Both experiments capture basic properties of the regional climate, like the southwest-northeast tilt in the North Atlantic storm track and the precipitation distribution across Europe. The westerlies and atmospheric moisture transport are stronger in RES than in CONTROL. Moreover, the tilt of the storm track in RES is reduced a little. These differences between the experiments result in increased precipitation in RES throughout most of the target domain, and especially over western Europe, to the south of the climatological maximum precipitation near the British Isles.
These differences appear physically plausible; we forego a more detailed analysis of the associated meteorological processes because of the comparative simplicity and low resolution of PlaSim. The results do illustrate how RES could be used in a more realistic model to investigate the physical processes driving extreme winter rainfall.

Conclusion
In this article we have demonstrated the applicability of RES algorithms for the study of extreme winter precipitation over Europe. We have applied a variant of the class of genealogical algorithms (Del Moral & Garnier, 2005). In our setup realizations in an ensemble simulation with sustained rainfall over Europe are cloned into multiple copies to favor the development of extreme winter rainfall.
As the results in Section 3 have shown, using the GPA RES algorithm described in Section 2.2.2 and Appendix B, we can sample much more efficiently realizations with extreme European winter precipitation. On the one hand, this allows to obtain more accurate estimates of rare event probabilities  7 of 11 at the same numerical cost. As demonstrated in Figure 4, for winter rainfall above approximately 0.2175 m, the RES becomes more efficient. On the other hand, estimates with the same accuracy can be obtained with less numerical effort. Our application demonstrates a combination of both, a relatively short RES of 61.45 years (compared to the 490 years of the control simulation) allows us to obtain probability estimates of rare events, with lower RE compared to the brute force estimates based on the control (see Figure 5).
Our results highlight the potential for further applications of RES in weather and climate science. Operational activities such as seasonal forecasting or event attribution are based on model experiments that are conceptually similar to the proof of concept presented in this study, albeit using much more sophisticated, and hence expensive, models. RES would permit these applications at lower cost and/or greater sophistication, hence allowing for a more in-depth investigation of the physical drivers of extremes. Moreover, RES permits to simulate events that are so rare/ extreme that they are simply not accessible within a brute-force approach while still being physically plausible within the chosen model and experiment. Extrapolations based on limit theorems, such as extreme value theory or large deviation theory, can provide probabilistic information, but not time-dependent geographical trajectories of individual realizations, as does GPA.
It should be stressed that the methodology used here is not specific to the general circulation model used. In the future, we plan to apply the method to a state-of-the-art, high-resolution coupled climate model. We also expect the method to work well for other geographical regions. Another possible avenue for future research include the testing of other types of weighting functions, taking into account the physical processes leading to rare events to sample them even more efficiently.  We exponentially tilt the distribution to ̃ with pdf̃( where  is a normalizing factor to ensure that the pdf integrates to 1.   . Notice that ̃ diverges for = 1 . The mean and standard deviation of ̃ change tõ As → 1 , both and 2 diverge.
If we wish to pick C so as to shift the mean to k standard deviations away from μ, we wish to have = + . Solving for C gives

Initiate
(d) If the model is not stochastic apply a random perturbation to make clones diverge over time. 3. Time-evolve ( ) −1 for i ∈ {1, …, N n−1 } under the model dynamics from t n−1 to t n = T 1 , resulting in ( ) ( ) with t ∈ (T 0 , T 1 In (Del Moral & Garnier, 2005), a central limit theorem result for the estimator result is derived for estimator B1. This result is important, because asymptotic confidence intervals can be derived from it. In the case of estimation of interval probabilities ℙ{ ∈ ( − Δ ∕2, + Δ ∕2)} with Δθ → 0, the leading order term in the variance can be estimated as̃(

Data Availability Statement
Jupyter notebooks and Julia scripts to run the rare event simulation described in this manuscript can be found at https://gitlab.act.reading.ac.uk/j.wouters/plasimgpa.jl. The Planet Simulator climate model can be downloaded from https://github.com/HartmutBorth/PLASIM (The version of PlaSim used in this study is Git commit ea4d2d8, dated 28/06/2017). This document has been drafted using GNU TEXMACS (van der Hoeven et al., N.D).