# Optimal Precursors Identification for North Atlantic Oscillation Using the Parallel Intelligence Algorithm

Authors: Bin Mu; Jing Li; Shijin Yuan; Xiaodan Luo; Guokun Dai
Journal: Scientific Programming
 (2022)

---

## Abstract

The North Atlantic Oscillation (NAO) with abnormal sea level pressure (SLP) differences influence the westerly wind strength and storm tracks in the North Atlantic, which further affects the winter climate in the northern hemisphere. The predictability of NAO has become an important area of climate research in recent years. The identification of optimal precursors (OPRs) would help to investigate the dynamics of atmospheric and oceanic motions, as well as the nonlinear characteristic. The conditional nonlinear optimal perturbation (CNOP) method has been widely used in research on the OPR of the climatic event to explore which kind of initial perturbation is most likely to trigger climate events. However, the previous works on NAO’s OPR are based on a simple ideal model, which cannot describe the evolution process of NAO. Moreover, the commonly used algorithms that rely on the adjoint model are not suitable for the large complicated numerical models that do not have corresponding adjoint models, like Community Earth System Model (CESM). To break through the limitation, this paper proposes the parallel principal component analysis (PCA)-based particle swarm optimization (PSO) and genetic algorithm (GA) hybrid algorithm (PGAPSO) algorithm to identify the OPR of NAO using CESM. For different initial conditions, OPRs identified by the proposed method can always make the basic state develop into NAO events with corresponding phases. Compared with other adjoint-free approaches, the proposed method has relatively better fitness values and more mature patterns. The results also reveal that the proposed method can avoid falling into a local optimum and has strong robustness. In terms of performance, the proposed method searches in feature space with a much lower number of dimensions, thus improving efficiency. In addition, the proposed method is accelerated with multiple parallel frameworks to enhance performance, and it achieves the speedup ratio of 40.2.


---

## Body

## 1. Introduction

The North Atlantic Oscillation (NAO) is the continuous fluctuation of the sea level pressure (SLP) field in the North Atlantic, which can be related to the interannual variation of the pressure over the Azores and Iceland. It represents large-scale alterations in the SLP differences between the subtropical and subpolar regions of the North Atlantic [1]. As the dominant variability mode of the atmospheric circulation in the northern hemisphere, the NAO is the result of complex nonlinear interactions between many spatiotemporal scales [2]. Previous research has shown that the NAO may be strictly linked to the anticyclonic (cyclonic) Rossby wave breaking [3] and can be viewed as a process with an e-folding time scale of about two weeks [4]. The NAO index (NAOI) is a quantitative indicator of the NAO, and its classical definition is the normalized difference between SLP over Iceland and over Azores [5]. The NAOI describes the intensity of the west wind drift in the North Atlantic and is closely related to the winter in Europe. The positive phase of the NAO (NAO+) is under a positive NAOI, while the negative phase of the NAO (NAO−) is under a negative NAOI. In the past decade, the fluctuation of the winter NAO has been quite extreme, and it has contributed significantly to the warm winter phenomenon throughout Europe, the cold weather in the Northwest Atlantic [6], the dipole precipitation pattern over northwest Europe and northwest Africa [7], and the surface temperature variation in the Iberian Peninsula and northeastern Europe [8], for instance.The NAO can be regarded as a nonlinear initial value problem [9], and its optimal precursor (OPR) is a kind of initial perturbation that is most likely to develop into NAO events [10]. The first kind of predictability problem about inaccurate initial conditions can be explored in the process of OPR's computation [11]. Since the initial condition has a significant influence on the predictability of the decadal variability [12], the OPR can help us to understand the dynamical processes of the NAO state transition. The physical mechanism for triggering the NAO event can be discovered by investigating the developing process of the OPR. Moreover, the sensitive areas determined by the spatial structure of the OPRs are beneficial to the intensive observations, thus improving the forecast accuracy of the NAO state transition. Furthermore, on the basis of the OPR, the target observation sensitive areas can also be detected. If the intensive observation is conducted in the target observation sensitive area, higher observation returns would be obtained. Hence, the research on OPR is of widespread scientific research value to study the physical mechanism and enhance forecast skills for the NAO. Further, exploring and optimizing the method for solving the OPR also has essential meanings. Although the phase and amplitude of NAO are affected by numerous factors [13–15], the characteristics of the NAO in the atmospherical process can be captured by the nonlinear models [16]. As a new generation of ocean-atmosphere coupled model, the Community Earth System Model (CESM) has presented the outstanding performance of NAO simulation, and the low-frequency variability of the NAO has also been well demonstrated [17, 18].The conditional nonlinear optimal perturbation (CNOP) is a mature method for solving OPRs. It describes the initial perturbation that causes the largest prediction error under a specific constraint condition at the prediction time. CNOP is appropriate for predictability studies of climate events with simulating nonlinear motions of oceans and atmospheres [19] and can be concluded as an optimization problem with constraints. CNOP approach was initially adopted to identify the OPRs of ENSO [20] and was gradually applied in research on the onset of blocking events [21], Kuroshio large meander [22], and Indian Ocean dipole events [23]. Recently, Jiang et al. [24] explored the OPRs that trigger the NAO events using CNOP, demonstrating that the amplitude induced by the self-interaction of perturbations in the onset of the NAO− is stronger than that in the onset of the NAO+. On this basis, Dai et al. [25] investigated the relationship between the OPR and optimally growing initial error (OGE) using CNOP. It was indicated that the two types of OGEs and the OPRs corresponding to the two types of NAO events have similar structures, and both of them can develop into dipole NAO anomaly patterns. These studies have proved that CNOP is a useful tool to investigate the onset of the NAO event. In their studies, the T21L3 quasigeostrophic global spectral model, which is a simple three-level model designed by Marshall and Molteni [26], was applied under ideal conditions. Due to the feature of the T21L3 model, they selected geopotential height as the characterized variable, and potential vorticity is the input variable. For solving CNOP, Jiang et al. [24] and Dai et al. [25] all used the spectral projected gradient 2 (SPG2) algorithm [27]. The SPG2 was designed to solve the minimum problem with restraints by determining the gradients of the cost function [28]. Several approaches with the same type have also been adopted to calculate the CNOP in the early years [29–31], such as the sequential quadratic programming (SQP) algorithm [32] and the limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm [33], etc. Since these algorithms rely on gradient information, the corresponding adjoint model needs to be called to obtain the gradients of the initial condition in the solving process.However, these adjoint-based methods would have a high probability to produce local CNOPs when the objective function has multiple extreme values and would fail with large initial disturbance or long prediction time due to the strong nonlinearity of the dynamical model [34]. Crucially, these traditional adjoint-based algorithms are not feasible to solve CNOP in complicated operational models that do not have an adjoint available [35]. In recent years, swarm intelligence algorithms have been gradually applied in the research of the CNOP [36, 37]. These algorithms determined the search direction from the position and fitness values of particles instead of gradients so that they can be extended to models without the adjoint model. It was also indicated that the swarm intelligence method still achieved global CNOP and had a shorter running time in the situation of larger initial perturbations, longer prediction times, multiple extrema values [34], and discontinuous objective functions [38]. Although these algorithms have shortened the runtime, it was still very time-consuming to search CNOP in the original dimensions. To enhance the performance, the researchers combined feature extraction strategies with intelligent algorithms, transforming the problems with high dimensions into the low-dimension space [39]. At present, the tentative application of intelligent algorithms based on feature extraction yielded concrete achievements. The principal component analysis based genetic algorithm (PCAGA) [40], the modified artificial bee colony algorithm (MABC) [41], the dynamic search fireworks algorithm with linearly decreased dimension number strategy (ld-dynFWA) [42], and PCA-based flower pollination (PCAFP) [43] have been successfully applied in the researches of tropical cyclone adaptive observations, El Niño-Southern Oscillation, and double-gyre variation, respectively. The CNOPs obtained by these methods had similar patterns and larger fitness values in comparison to the adjoint method. It was indicated that the PCA-based intelligence algorithm is appropriate for solving CNOP in high-dimensional numerical models, especially the models without the adjoint model, like CESM.The objective of this paper is to find the OPRs which can produce the NAO anomaly pattern and explore the effect of the nonlinear process. We attempt to answer the following questions. (1) How to solve the CNOP and identify the OPR effectively and efficiently? (2) What is the intensity of maximum uncertainty of NAO prediction triggered by the OPRs? (3) Which kind of spatial structure would cause these two phases of NAO events? (4) How does the SLP field change during the prediction period after superimposing OPRs? The topic is studied using the CESM, which is an ocean-atmosphere coupled model without an adjoint model. For the first problem, the PCA-based intelligence algorithm is adopted to identify the OPR using CESM. In previous studies related to CNOP, it was found that particle swarm optimization (PSO) has a significant advantage but still suffers from the disadvantage of being trapped in a local optimum, which can be effectively mitigated by coevolutionary strategies with other intelligence algorithms [44, 45]. Thus, the hybrid algorithm named PCA-based GA and PSO (PGAPSO) is proposed to solve the CNOP. To further improve the efficiency of the algorithm, the parallel scheme using the Message Passing Interface (MPI) and the Compute Unified Device Architecture (CUDA) is utilized to enhance the time performance of both the algorithm and CESM.The structure of this paper is organized as follows: Section2 describes the CESM, the CNOP method, the PGAPSO algorithm, and the parallelization technique. Experiments and results are displayed in Section 3. This paper ends with a conclusion and future work in Section 4.


## 2. Materials and Methods

### 2.1. CESM

The CESM [46] is a new generation of fully coupled climate model developed in 2010. It has been widely used to simulate the carbon cycle [47], ocean currents [48], soil moisture [49], precipitation [50], and other climate phenomena. CESM is composed of seven geophysical models, respectively, Atmosphere (ATM), Sea-ice (ICE), Land (LND), River-runoff (ROF), Ocean (OCN), Land-ice (GLC), and Ocean-wave (WAV, stub only). The CESM also has a coupler (CPL) that coordinates the time evolution of geophysical models and delivers information between these components.CESM is currently updated to version 2.2, and the cloud layers unified by binormals (CLUBB) scheme replace the parameterization scheme in the previous stable version 1.2.2. The execution time is increased correspondingly with these modifications. It even needs to take several hours to simulate a NAO process during the few days, which is much longer than that of CESM1. To test their simulation capability of the NAO, CESM 2.2.0 and CESM 1.2.2 are both executed for a 90-day simulation in winter starting from the same reference time. The correlation coefficient of the 90-day NAOI between CESM 1.2.2 and CESM 2.2.0 is 0.83, which demonstrates a similar simulation.Because of the reasonable running time and simulation capability for NAO simulation, the atmospherical component in CESM 1.2.2 is selected in this work. The Community Atmosphere Model (CAM) version 5.3, which is a global atmospherical general circulation model developed from the NCAR CCM3, is released as the atmosphere component of CESM 1.2. The CAM incorporates an interactive aerosol model where aerosols interact with the tropospheric chemistry. The Community Land Model (CLM) is adopted as the active land model in CESM, and the River Transport Model (RTM) is used to route total runoff from the land model to oceans. The component setF_2000 that includes an active atmospherical component (CAM), active land component (CLM), active river-runoff (RTM) component, and sea-ice component (the sea-ice component, CICE) in the prescribed mode is selected. The GLC and WAV models are in stub mode (SGLC and SWAV), which only satisfy interface requirements and are not needed for configuration. F_2000 also includes the data ocean component (DOCN). It reads SSTs from input datasets and can run as a pure data model. The models in component set F_2000 are shown in Figure 1.Figure 1
The models in component setF_2000; active models are in green, data models are in blue, and stub models are in orange.The NAO cases are simulated in the winter (DJF), and the parameternhtfrq is set to −24, which denotes the daily average. The experiments are performed on a resolution of f09_g16 with an approximate grid spacing of 0.9°×1.25°, including 26 levels in the vertical. Perturbations are superimposed on the basic state of the Arctic region (north of 60°N, see Figure 2) and contain six variables listed in Table 1. The size of U, V, T, and Q is 26 (layer) × 32 (latitude) × 288 (longitude), and the size of Π and Φ is 32 × 288. The NAOI increment is selected as the fitness value and is calculated by the output variable SLP. The region of SLP is mainly located between 20°N and 80°N and between 90°W and 40°E, as shown in Figure 3. The vector of SLP only has one layer and consists of 65 × 105 grids.Figure 2
The regions of perturbations and the North Atlantic sector where the NAO events mainly occur. The perturbations superimpose on the Arctic section between60°N and 90°N (shown in red). The region of NAO is between 20°N and 80°N and between 90°W and 40°E (shown in yellow), and it is also the scope of the SLP variable.Table 1
The related variables included in the perturbations.
Variable nameDescriptionUnitsUZonal windm/sVMeridional windm/sTTemperatureKQSpecific humiditykg/kgΠSurface pressurePaΦSurface geopotentialm2/s2Figure 3
The first mode of the EOF with SLP anomaly field concentrated in the North Atlantic region between90°W–40°E and 20°N–80°N.


### 2.2. CNOP

The CNOP is a natural extension of the linear singular vector into the nonlinear regime and is proposed to study predictability problems of weather and climate in numerical models [51, 52]. The CNOP can represent the initial perturbation that can trigger the largest uncertainty at the prediction time. The OPR of NAO is a kind of CNOP that makes the basic state evolve into the NAO events. Suppose the nonlinear model can be briefly described as follows:(1)∂S∂t+FS=0,S|t=t0=S0,where S denotes the state vector and S0 is the basic state at the initial time t0. F is the nonlinear partial differential operator of the model. Equation (1) has the following discrete form:(2)St=Mt0⟶tS0,where Mt0⟶t represents the nonlinear propagator that “propagates” the initial state from time t0 to the prediction time t. St is the state vector at time t. If an initial perturbation s0 on the basic state is superimposed, the development of the initial perturbation would be(3)ΔS=St′−St=Mt0⟶tS0+s0−Mt0⟶tS0,where St′ denotes the final state with superposing perturbations, while ΔS stands for the development of the final state, which is the difference between perturbation state and reference state. The OPR refers to the initial perturbation that can make the objective function achieve the maximum (minimum) under the constraint condition at prediction time. In this work, the NAOI difference between the perturbation state and the reference state is selected as the objective function J:(4)Js0∗NAO+=maxcs0≤ϵJs0=NAOINAO+CNOP−NAOIrefer,(5)Js0∗NAO−=mincs0≤ϵJs0=NAOINAO−CNOP−NAOIrefer,where s0 is the initial perturbation and consists of physical variables listed in Table 1. According to the definition of the OPR, the perturbation s0∗ (NAO+) is the OPR of the NAO+ and makes J achieve the maximum, whereas s0∗ (NAO−) is the OPR of the NAO− and makes J achieve the minimum. NAOINAO+CNOP refers to the highest positive NAOI, which is triggered by CNOP-typed perturbation, while NAOINAO−CNOP is the lowest negative NAOI caused by CNOP. These two types of CNOPs are solved with opposite search directions. cs0 denotes the constraint function of perturbation s0, and ϵ denotes the constraint condition of the OPRs, ensuring the perturbation within a reasonable range. The constraint condition is consulted from the related works of the sensitive area identification for tropical cyclone [40] owing to the same variables and is set to 10% of the summation of the kinetic energy of basic state in the vertical coordinate σ and verification areas D (i.e., the area north of 60°N):(6)cs0=1D∫D∫01u2+v2+CpTrpt2+RaTrππr2dσdDϵ=10%∗1D∫D∫01U02+V02+CpTrT02+RaTrΠ0πr2dσdD,where u, v, pt, and π are perturbation vectors of zonal wind, meridional wind, temperature, and surface pressure, respectively, and U0, V0, T0, and Π0 are the initial conditions of these variables. For all six variables in perturbation s0, u, v, pt, and π are constrained via equation (5), and the other two variables q and ϕ are constrained by the budget check and balance check built into CESM. Cp is the specific heat at the constant pressure which is set to 1005.7 J⋅kg−1⋅K−1 and Tr is the reference temperature with a value of 270K. Ra denotes the ideal gas constant, and its value is set to 287.05 J⋅kg−1⋅K−1. πr is the static reference pressure with a value of 1000 hPa. The constraint condition is static for a certain case since the constraint should specify the range of initial perturbations before the simulation starts. Thus, the constraint can only be defined by physical quantities in the initial field and is a solid hypersphere with radius ϵ. It is necessary to ensure u, v, pt, and π are all located inside the hypersphere or on the surface of the hypersphere to make the model input reasonable. If the initial perturbation is out of the scope of constraint, the projection operation is required to map the perturbation inside the proper range:(7)s=s∗ϵd,d=s12+s22+⋯+sn2,where s is the initial perturbation and n denotes the dimension number of the solution space.NAOI can be defined in several ways, like normalized SLP difference between35°N and 65°N, and so on. In this experiment, a blocking indicator proposed by Liu [53] is chosen to quantify the extent of the NAO events. The NAOI is defined as the surface conserving projection of the SLP field on the NAO anomaly pattern:(8)NAOI=SLPNAO,SLPdSLPNAO,SLPNAO,specifically, NAOI=∑g=1lat∑h=1lonSLPNAOgh∗coslatg∗π/180∑g=1lat∑h=1lonSLPNAOgh∗SLPNAOgh∗ψ,where SLPd is obtained by subtracting the climatological mean from SLP output and ⋅,⋅ denotes the inner product operation of vectors. The details for the inner product are shown in the following line. lat and lon denote the region of the NAO pattern, and lat=65, lon=105. latg stands for the latitude of the current point, and the cosine-related coefficient is to ensure the homogeneity of the grid. ψ is the normalization coefficient and is set to 31914.543 to make the NAOI obey the normal distribution. SLPNAO denotes the NAO anomaly pattern decomposed by the empirical orthogonal function (EOF) analysis, which is illustrated in Figure 3. The NAO spatial pattern is manifested as a typical meridional dipole mode, including the Iceland low pressure along with the North Atlantic subtropical high. In Figure 3, it is a pattern of the NAO+, presenting the mode with the negative anomalies in high latitude and the positive anomalies in low latitude.Combining equations (4), (5), and (8), the objective function is described as follows:(9)Js0=ΔNAOI=Mt0⟶tS0+s0−Mt0⟶tS0,SLPNAOSLPNAO,SLPNAO.In summary, the objective function is the projection of the SLP field difference between the final state and the reference state on the NAO anomaly pattern, and the OPR ofNAO+ (NAO−) can be obtained by solving the maximum (minimum) of the objective function J. In general, NAOI is in the range of [−3, 3]. To simplify the problem, assume that an NAO event occurs when Δ NAOI > 1.0 (NAO+ event) or Δ NAOI < −1.0 (NAO− event).


### 2.3. PGAPSO

Under the resolution of0.9°×1.25°, the total dimensions of variables involved in the objective function are 4 (U, V, T, and Q) × 26 (layer) × 32 (latitude) × 288 (longitude) + 2 (Π, Φ) × 32 (latitude) × 288 (longitude) = 976896. It is difficult for the algorithm to search for the optimal solution in such high dimensions. Thus, the scale of the solution domain is reduced via eigenvector decomposition using the PCA method.First, 900 original samples are obtained after running 10-year integration (only in winter) on a daily average using CESM. Owing to their different physical units, Z-Score normalization is needed to be performed for each variable:(10)Pv=Pv−Pv¯σvv=1,2,…,c,where c refers to the number of variables; that is, c=6. Pv¯ is the mean value of the variable v, and the σv denotes the standard deviation.After that, each piece of sample is reshaped into a vector with one dimension (1× 976896), these samples constitute a matrix with a size of 900 × 976896. Then, each sample is subtracted from the mean values of these samples for centering, and the obtained sample is weighted according to the area of the grid:(11)Pij′=Pij−1n∑i=1nPij∗coslatji=1,2,…,n;j=0,1,…,r,where Pij denotes the jth element of the ith sample data, n is the number of the samples, and r is the length for each sample. n equals 900, and r equals 976896. latj is the corresponding latitude of the jth element, and the area weight is calculated approximately via the cosine value of the latitude. Then, the eigenvalues (λ1, …, λn) and eigenvectors of the covariance matrix PPT are calculated to obtain principal components:(12)PPTL=LΣ,where L is the eigenvector matrix and Σ is a diagonal matrix whose entries in the main diagonal are the corresponding eigenvalues. The top m columns of the eigenvectors L sorted by their eigenvalues are selected as the principal components C with a size of m×n. The variance contribution rate λi/∑k=1nλk reflects the amount of information that the component Ci contains, and the value of m is determined by the accumulative variance ratio ∑k=1mλk/∑k=1nλk.To obtain the extremum of the objective function, a hybrid algorithm improved from two efficient algorithms named PSO [54] and GA is adopted. Since PSO has a relatively better performance in solving CNOP [44, 45], it is adopted as the main body of the PGAPSO. The PSO conducts search via particles, and each particle contains two properties: position (X) and speed (V). The position is a group of solutions and is updated with the help of speed. The position with the optimal fitness value is the OPR that we desire. The flow of the PSO can be briefly described as follows:(1)
Initialize the speed and position of each particle with random values.(2)
For each particlef, the position vector is in reduced space, so the position vector needs to be restored into original space via Xf′=Xf⋅C. Then, superpose the perturbation Xf′ on the basic state. When the model integration is finished, calculate the fitness value of each particle through equation (9) and record its optimal position (Xpbf) along with the optimal global position (Xgb). Xpbf denotes the position of particle f with its optimal fitness value, while Xgb stands for the optimal position of all particles.(3)
Update the position and speed of each particle. The position with the best fitness value is searched by tracing both the individual optimal position and the optimal global position:(13)Vfk+1=ωkVfk+c1r1Xpbfk−Xfk+c2r2Xgbk−Xfk,Xfk+1=Xfk+Vfk+1f=1,…,N,whereVfk+1 is the speed of particle f for step k+1 and Xpbfk refers to the best position of particle f up to step k. Xpbk represents the best position of the entire swarm in k steps. N is the total number of particles. c1 and c2 are the self-awareness coefficient and social-awareness coefficient, respectively. r1 and r2 are random parameters.ωk is the inertia weight parameter. To improve the convergence speed, the algorithm is designed to transform from the global search to the local search as the number of iterations increases. After several iterations, the particle speed is limited via the linearly decreasing ωk:(14)ωk=ωmax−ωmax−ωminitermax∗k.(4)
Repeat steps as mentioned above until satisfying the terminal condition.Based on the strategy of PSO, GA further optimizes the search process. To avoid trapping in local optimum, the particles are fed into GA for further search when the trend of convergence is evident in the PSO. The mean speed normVk2¯ for step k is adopted to measure the convergence of the PSO:(15)Vk2¯=∑f=1NVfk2N,whereVfk2=∑i=1nvi2.Thus, when the mean speed norm is smaller than the fusion thresholdξ, the particles are updated by the GA procedure.As a metaheuristic algorithm, the GA derives from natural selection [55]. The main operations of GA include selection, crossover, and mutation. Particles are selected according to their fitness value to breed a new generation. The selection probability for each individual is equal to the ratio of its fitness value to the sum of the fitness value population:(16)ps=JsXf′∑JsX′.Obviously, individuals with larger fitness values have a greater probability of being selected. Furthermore, the two selected individualsa and b with positions of Xax0,…,xs,…,xe,…,xm and Xbx0,…,xs,…,xe,…,xm have a pc probability of cross-recombining to generate two recombinant individuals with position of Xa′ and Xb′:(17)Xa′xs,…,xe=Xbxs,…,xe,Xb′xs,…,xe=Xaxs,…,xe,xs and xe denote the starting point and ending point of the crossover, respectively, satisfying 0≤s<e<m. a and b exchange their position vector from xs to xe. After the crossover, two new individuals a′ and b′ replace the original individuals for the calculation of fitness values.To increase randomness and improve its global search ability, individuals have apm probability of mutation. Here, the variation is limited to a single-point mutation, where a random point within [0, m − 1] is replaced by a random number to form a new individual. The fitness value of each new particle is compared against its parents, and the best position is recorded. If the new individual has a better fitness value than the global best position, the global best position (Xgbk) would be replaced. Then, the speed and position of particles are updated using equation (13), and the final global fitness value is obtained until the iter reaches itermax.The parameter settings of the PGAPSO are listed in Table2.Table 2
The definition and value of controls parameter in the PGAPSO.
Parameter nameDescriptionValuenParticle number30mDimension of feature space50itermaxMaximum number of iterations50ωmaxUpper limit of inertia weight parameter1ωminLower limit of inertia weight parameter0.5c1Self-awareness coefficient2c2Social-awareness coefficient2r1Random parameterRandom value in [0, 1]r2Random parameterRandom value in [0, 1]ξFusion threshold20pcCrossover probability0.3pmMutation probability0.1


### 2.4. Parallelization

The computation of CNOP in CESM is quite time-consuming. With 48 CPU cores, 30 particles, and 50 iterations, it takes about 8.33 days to obtain the OPRs in the serial program. Multiple parallel techniques and frameworks are adopted to improve the performance of this work.The role of the CAM component in CESM is to simulate the variation of atmosphere, and the integration process involves plenty of matrix operations with high-dimensional input data. Thus, the model simulation becomes the bottleneck of the program performance. Recently, the Graphics Processing Unit (GPU) has been widely used in accelerating numerical models. Since GPU is suitable for large-scale parallel computing, it can significantly improve the execution performance of numerical calculations in climate models.Previous work results have shown that the performance of the optimized subroutine can be improved by the GPU technique substantially [56–59]. In this work, subroutine radclwmx and radabs both execute longer according to the performance evaluation using perf. The subroutine radclwmx adopts the broadband absorption method to calculate the longwave emissivity and boundary fluxes under the assumption that random overlapping clouds have variable cloud emissivity. The subroutine radabs computes the absorptivities for gases like H2O, CO2, O3, CH4, N2O, CCl3F, and CCl2F2. These two subroutines are called with high frequency and are rewritten with the CUDA interface. Simultaneously, kernel directives and OpenACC directives are adopted to simplify specific operations on the device, and multiple asynchronous streams are also created to overlap function execution and data replication. For the compilation phase, the option -O4 is selected to perform the optimization at the highest level. The option -fast and -fastsse are also utilized to launch the 64-bit Single Instruction Multiple Data (SIMD) instruction and implement cache alignment and flush.In each iteration of PGAPSO, the calculation of the fitness value for each particle is relatively independent. According to its features, multiprocess techniques are suitable for executing these tasks concurrently. Here, MPI is adopted as the parallel framework to accelerate the algorithm. The pseudocode describes the flow of the parallel PGAPSO in Table3.Table 3
The pseudocode of the procedures in PGAPSO.
Algorithm 1 pseudocode for PGAPSO1: transform training data through PCA to obtain principal componentsLd with m dimensions2: initialize population3: FORiter = 1 to iter = itermax DO4:  FOR each processi DO5:   restore solution matrixXi into original space via Xi∗Ld6:   calculate fitness valueJs0=FMt0⟶tS0+P∗Ld under the constraint ϵ7:  END FOR8:  IF norm of particle speed≤ξ THEN9:   select individuals according tops=JsXi′/∑JsX′10:   WHILE count≤ N/2 DO11:    crossover with probabilitypc12:   END WHILE13:   WHILE count≤ N DO14:    mutate with probabilitypm15:   END WHILE16:   compare fitness values between new generation and parent individuals17:  END IF18:  update particle speed viaXpbk and Xgbk19:  update particle position20: END FORThese processes are divided into two groups: one is a master process; the others are slave processes. Figure4 illustrates the parallel architecture of PGAPSO for solving CNOP. At each iteration, the master process allocates calculation tasks to slave processes. For each slave process, perturbation under constraint condition is generated and is superimposed on the basic state of CESM. CESM, which is paralleled with MPI, OpenMP, and CUDA, is called to perform the integration. The fitness values are calculated via equation (9) and are gathered by the master process. After that, the master process broadcasts the optimal global value to slave processes via MPI. If particles’ speed norm is less than the threshold value, the crossover and mutation operations are conducted. Then, the master process compares the current particles and updates the information of optimal particles at the end of each iteration. With the help of MPI, the performance of PGAPSO can be significantly enhanced.Figure 4
The parallel architecture of PGAPSO for solving CNOP.


## 2.1. CESM

The CESM [46] is a new generation of fully coupled climate model developed in 2010. It has been widely used to simulate the carbon cycle [47], ocean currents [48], soil moisture [49], precipitation [50], and other climate phenomena. CESM is composed of seven geophysical models, respectively, Atmosphere (ATM), Sea-ice (ICE), Land (LND), River-runoff (ROF), Ocean (OCN), Land-ice (GLC), and Ocean-wave (WAV, stub only). The CESM also has a coupler (CPL) that coordinates the time evolution of geophysical models and delivers information between these components.CESM is currently updated to version 2.2, and the cloud layers unified by binormals (CLUBB) scheme replace the parameterization scheme in the previous stable version 1.2.2. The execution time is increased correspondingly with these modifications. It even needs to take several hours to simulate a NAO process during the few days, which is much longer than that of CESM1. To test their simulation capability of the NAO, CESM 2.2.0 and CESM 1.2.2 are both executed for a 90-day simulation in winter starting from the same reference time. The correlation coefficient of the 90-day NAOI between CESM 1.2.2 and CESM 2.2.0 is 0.83, which demonstrates a similar simulation.Because of the reasonable running time and simulation capability for NAO simulation, the atmospherical component in CESM 1.2.2 is selected in this work. The Community Atmosphere Model (CAM) version 5.3, which is a global atmospherical general circulation model developed from the NCAR CCM3, is released as the atmosphere component of CESM 1.2. The CAM incorporates an interactive aerosol model where aerosols interact with the tropospheric chemistry. The Community Land Model (CLM) is adopted as the active land model in CESM, and the River Transport Model (RTM) is used to route total runoff from the land model to oceans. The component setF_2000 that includes an active atmospherical component (CAM), active land component (CLM), active river-runoff (RTM) component, and sea-ice component (the sea-ice component, CICE) in the prescribed mode is selected. The GLC and WAV models are in stub mode (SGLC and SWAV), which only satisfy interface requirements and are not needed for configuration. F_2000 also includes the data ocean component (DOCN). It reads SSTs from input datasets and can run as a pure data model. The models in component set F_2000 are shown in Figure 1.Figure 1
The models in component setF_2000; active models are in green, data models are in blue, and stub models are in orange.The NAO cases are simulated in the winter (DJF), and the parameternhtfrq is set to −24, which denotes the daily average. The experiments are performed on a resolution of f09_g16 with an approximate grid spacing of 0.9°×1.25°, including 26 levels in the vertical. Perturbations are superimposed on the basic state of the Arctic region (north of 60°N, see Figure 2) and contain six variables listed in Table 1. The size of U, V, T, and Q is 26 (layer) × 32 (latitude) × 288 (longitude), and the size of Π and Φ is 32 × 288. The NAOI increment is selected as the fitness value and is calculated by the output variable SLP. The region of SLP is mainly located between 20°N and 80°N and between 90°W and 40°E, as shown in Figure 3. The vector of SLP only has one layer and consists of 65 × 105 grids.Figure 2
The regions of perturbations and the North Atlantic sector where the NAO events mainly occur. The perturbations superimpose on the Arctic section between60°N and 90°N (shown in red). The region of NAO is between 20°N and 80°N and between 90°W and 40°E (shown in yellow), and it is also the scope of the SLP variable.Table 1
The related variables included in the perturbations.
Variable nameDescriptionUnitsUZonal windm/sVMeridional windm/sTTemperatureKQSpecific humiditykg/kgΠSurface pressurePaΦSurface geopotentialm2/s2Figure 3
The first mode of the EOF with SLP anomaly field concentrated in the North Atlantic region between90°W–40°E and 20°N–80°N.


## 2.2. CNOP

The CNOP is a natural extension of the linear singular vector into the nonlinear regime and is proposed to study predictability problems of weather and climate in numerical models [51, 52]. The CNOP can represent the initial perturbation that can trigger the largest uncertainty at the prediction time. The OPR of NAO is a kind of CNOP that makes the basic state evolve into the NAO events. Suppose the nonlinear model can be briefly described as follows:(1)∂S∂t+FS=0,S|t=t0=S0,where S denotes the state vector and S0 is the basic state at the initial time t0. F is the nonlinear partial differential operator of the model. Equation (1) has the following discrete form:(2)St=Mt0⟶tS0,where Mt0⟶t represents the nonlinear propagator that “propagates” the initial state from time t0 to the prediction time t. St is the state vector at time t. If an initial perturbation s0 on the basic state is superimposed, the development of the initial perturbation would be(3)ΔS=St′−St=Mt0⟶tS0+s0−Mt0⟶tS0,where St′ denotes the final state with superposing perturbations, while ΔS stands for the development of the final state, which is the difference between perturbation state and reference state. The OPR refers to the initial perturbation that can make the objective function achieve the maximum (minimum) under the constraint condition at prediction time. In this work, the NAOI difference between the perturbation state and the reference state is selected as the objective function J:(4)Js0∗NAO+=maxcs0≤ϵJs0=NAOINAO+CNOP−NAOIrefer,(5)Js0∗NAO−=mincs0≤ϵJs0=NAOINAO−CNOP−NAOIrefer,where s0 is the initial perturbation and consists of physical variables listed in Table 1. According to the definition of the OPR, the perturbation s0∗ (NAO+) is the OPR of the NAO+ and makes J achieve the maximum, whereas s0∗ (NAO−) is the OPR of the NAO− and makes J achieve the minimum. NAOINAO+CNOP refers to the highest positive NAOI, which is triggered by CNOP-typed perturbation, while NAOINAO−CNOP is the lowest negative NAOI caused by CNOP. These two types of CNOPs are solved with opposite search directions. cs0 denotes the constraint function of perturbation s0, and ϵ denotes the constraint condition of the OPRs, ensuring the perturbation within a reasonable range. The constraint condition is consulted from the related works of the sensitive area identification for tropical cyclone [40] owing to the same variables and is set to 10% of the summation of the kinetic energy of basic state in the vertical coordinate σ and verification areas D (i.e., the area north of 60°N):(6)cs0=1D∫D∫01u2+v2+CpTrpt2+RaTrππr2dσdDϵ=10%∗1D∫D∫01U02+V02+CpTrT02+RaTrΠ0πr2dσdD,where u, v, pt, and π are perturbation vectors of zonal wind, meridional wind, temperature, and surface pressure, respectively, and U0, V0, T0, and Π0 are the initial conditions of these variables. For all six variables in perturbation s0, u, v, pt, and π are constrained via equation (5), and the other two variables q and ϕ are constrained by the budget check and balance check built into CESM. Cp is the specific heat at the constant pressure which is set to 1005.7 J⋅kg−1⋅K−1 and Tr is the reference temperature with a value of 270K. Ra denotes the ideal gas constant, and its value is set to 287.05 J⋅kg−1⋅K−1. πr is the static reference pressure with a value of 1000 hPa. The constraint condition is static for a certain case since the constraint should specify the range of initial perturbations before the simulation starts. Thus, the constraint can only be defined by physical quantities in the initial field and is a solid hypersphere with radius ϵ. It is necessary to ensure u, v, pt, and π are all located inside the hypersphere or on the surface of the hypersphere to make the model input reasonable. If the initial perturbation is out of the scope of constraint, the projection operation is required to map the perturbation inside the proper range:(7)s=s∗ϵd,d=s12+s22+⋯+sn2,where s is the initial perturbation and n denotes the dimension number of the solution space.NAOI can be defined in several ways, like normalized SLP difference between35°N and 65°N, and so on. In this experiment, a blocking indicator proposed by Liu [53] is chosen to quantify the extent of the NAO events. The NAOI is defined as the surface conserving projection of the SLP field on the NAO anomaly pattern:(8)NAOI=SLPNAO,SLPdSLPNAO,SLPNAO,specifically, NAOI=∑g=1lat∑h=1lonSLPNAOgh∗coslatg∗π/180∑g=1lat∑h=1lonSLPNAOgh∗SLPNAOgh∗ψ,where SLPd is obtained by subtracting the climatological mean from SLP output and ⋅,⋅ denotes the inner product operation of vectors. The details for the inner product are shown in the following line. lat and lon denote the region of the NAO pattern, and lat=65, lon=105. latg stands for the latitude of the current point, and the cosine-related coefficient is to ensure the homogeneity of the grid. ψ is the normalization coefficient and is set to 31914.543 to make the NAOI obey the normal distribution. SLPNAO denotes the NAO anomaly pattern decomposed by the empirical orthogonal function (EOF) analysis, which is illustrated in Figure 3. The NAO spatial pattern is manifested as a typical meridional dipole mode, including the Iceland low pressure along with the North Atlantic subtropical high. In Figure 3, it is a pattern of the NAO+, presenting the mode with the negative anomalies in high latitude and the positive anomalies in low latitude.Combining equations (4), (5), and (8), the objective function is described as follows:(9)Js0=ΔNAOI=Mt0⟶tS0+s0−Mt0⟶tS0,SLPNAOSLPNAO,SLPNAO.In summary, the objective function is the projection of the SLP field difference between the final state and the reference state on the NAO anomaly pattern, and the OPR ofNAO+ (NAO−) can be obtained by solving the maximum (minimum) of the objective function J. In general, NAOI is in the range of [−3, 3]. To simplify the problem, assume that an NAO event occurs when Δ NAOI > 1.0 (NAO+ event) or Δ NAOI < −1.0 (NAO− event).


## 2.3. PGAPSO

Under the resolution of0.9°×1.25°, the total dimensions of variables involved in the objective function are 4 (U, V, T, and Q) × 26 (layer) × 32 (latitude) × 288 (longitude) + 2 (Π, Φ) × 32 (latitude) × 288 (longitude) = 976896. It is difficult for the algorithm to search for the optimal solution in such high dimensions. Thus, the scale of the solution domain is reduced via eigenvector decomposition using the PCA method.First, 900 original samples are obtained after running 10-year integration (only in winter) on a daily average using CESM. Owing to their different physical units, Z-Score normalization is needed to be performed for each variable:(10)Pv=Pv−Pv¯σvv=1,2,…,c,where c refers to the number of variables; that is, c=6. Pv¯ is the mean value of the variable v, and the σv denotes the standard deviation.After that, each piece of sample is reshaped into a vector with one dimension (1× 976896), these samples constitute a matrix with a size of 900 × 976896. Then, each sample is subtracted from the mean values of these samples for centering, and the obtained sample is weighted according to the area of the grid:(11)Pij′=Pij−1n∑i=1nPij∗coslatji=1,2,…,n;j=0,1,…,r,where Pij denotes the jth element of the ith sample data, n is the number of the samples, and r is the length for each sample. n equals 900, and r equals 976896. latj is the corresponding latitude of the jth element, and the area weight is calculated approximately via the cosine value of the latitude. Then, the eigenvalues (λ1, …, λn) and eigenvectors of the covariance matrix PPT are calculated to obtain principal components:(12)PPTL=LΣ,where L is the eigenvector matrix and Σ is a diagonal matrix whose entries in the main diagonal are the corresponding eigenvalues. The top m columns of the eigenvectors L sorted by their eigenvalues are selected as the principal components C with a size of m×n. The variance contribution rate λi/∑k=1nλk reflects the amount of information that the component Ci contains, and the value of m is determined by the accumulative variance ratio ∑k=1mλk/∑k=1nλk.To obtain the extremum of the objective function, a hybrid algorithm improved from two efficient algorithms named PSO [54] and GA is adopted. Since PSO has a relatively better performance in solving CNOP [44, 45], it is adopted as the main body of the PGAPSO. The PSO conducts search via particles, and each particle contains two properties: position (X) and speed (V). The position is a group of solutions and is updated with the help of speed. The position with the optimal fitness value is the OPR that we desire. The flow of the PSO can be briefly described as follows:(1)
Initialize the speed and position of each particle with random values.(2)
For each particlef, the position vector is in reduced space, so the position vector needs to be restored into original space via Xf′=Xf⋅C. Then, superpose the perturbation Xf′ on the basic state. When the model integration is finished, calculate the fitness value of each particle through equation (9) and record its optimal position (Xpbf) along with the optimal global position (Xgb). Xpbf denotes the position of particle f with its optimal fitness value, while Xgb stands for the optimal position of all particles.(3)
Update the position and speed of each particle. The position with the best fitness value is searched by tracing both the individual optimal position and the optimal global position:(13)Vfk+1=ωkVfk+c1r1Xpbfk−Xfk+c2r2Xgbk−Xfk,Xfk+1=Xfk+Vfk+1f=1,…,N,whereVfk+1 is the speed of particle f for step k+1 and Xpbfk refers to the best position of particle f up to step k. Xpbk represents the best position of the entire swarm in k steps. N is the total number of particles. c1 and c2 are the self-awareness coefficient and social-awareness coefficient, respectively. r1 and r2 are random parameters.ωk is the inertia weight parameter. To improve the convergence speed, the algorithm is designed to transform from the global search to the local search as the number of iterations increases. After several iterations, the particle speed is limited via the linearly decreasing ωk:(14)ωk=ωmax−ωmax−ωminitermax∗k.(4)
Repeat steps as mentioned above until satisfying the terminal condition.Based on the strategy of PSO, GA further optimizes the search process. To avoid trapping in local optimum, the particles are fed into GA for further search when the trend of convergence is evident in the PSO. The mean speed normVk2¯ for step k is adopted to measure the convergence of the PSO:(15)Vk2¯=∑f=1NVfk2N,whereVfk2=∑i=1nvi2.Thus, when the mean speed norm is smaller than the fusion thresholdξ, the particles are updated by the GA procedure.As a metaheuristic algorithm, the GA derives from natural selection [55]. The main operations of GA include selection, crossover, and mutation. Particles are selected according to their fitness value to breed a new generation. The selection probability for each individual is equal to the ratio of its fitness value to the sum of the fitness value population:(16)ps=JsXf′∑JsX′.Obviously, individuals with larger fitness values have a greater probability of being selected. Furthermore, the two selected individualsa and b with positions of Xax0,…,xs,…,xe,…,xm and Xbx0,…,xs,…,xe,…,xm have a pc probability of cross-recombining to generate two recombinant individuals with position of Xa′ and Xb′:(17)Xa′xs,…,xe=Xbxs,…,xe,Xb′xs,…,xe=Xaxs,…,xe,xs and xe denote the starting point and ending point of the crossover, respectively, satisfying 0≤s<e<m. a and b exchange their position vector from xs to xe. After the crossover, two new individuals a′ and b′ replace the original individuals for the calculation of fitness values.To increase randomness and improve its global search ability, individuals have apm probability of mutation. Here, the variation is limited to a single-point mutation, where a random point within [0, m − 1] is replaced by a random number to form a new individual. The fitness value of each new particle is compared against its parents, and the best position is recorded. If the new individual has a better fitness value than the global best position, the global best position (Xgbk) would be replaced. Then, the speed and position of particles are updated using equation (13), and the final global fitness value is obtained until the iter reaches itermax.The parameter settings of the PGAPSO are listed in Table2.Table 2
The definition and value of controls parameter in the PGAPSO.
Parameter nameDescriptionValuenParticle number30mDimension of feature space50itermaxMaximum number of iterations50ωmaxUpper limit of inertia weight parameter1ωminLower limit of inertia weight parameter0.5c1Self-awareness coefficient2c2Social-awareness coefficient2r1Random parameterRandom value in [0, 1]r2Random parameterRandom value in [0, 1]ξFusion threshold20pcCrossover probability0.3pmMutation probability0.1


## 2.4. Parallelization

The computation of CNOP in CESM is quite time-consuming. With 48 CPU cores, 30 particles, and 50 iterations, it takes about 8.33 days to obtain the OPRs in the serial program. Multiple parallel techniques and frameworks are adopted to improve the performance of this work.The role of the CAM component in CESM is to simulate the variation of atmosphere, and the integration process involves plenty of matrix operations with high-dimensional input data. Thus, the model simulation becomes the bottleneck of the program performance. Recently, the Graphics Processing Unit (GPU) has been widely used in accelerating numerical models. Since GPU is suitable for large-scale parallel computing, it can significantly improve the execution performance of numerical calculations in climate models.Previous work results have shown that the performance of the optimized subroutine can be improved by the GPU technique substantially [56–59]. In this work, subroutine radclwmx and radabs both execute longer according to the performance evaluation using perf. The subroutine radclwmx adopts the broadband absorption method to calculate the longwave emissivity and boundary fluxes under the assumption that random overlapping clouds have variable cloud emissivity. The subroutine radabs computes the absorptivities for gases like H2O, CO2, O3, CH4, N2O, CCl3F, and CCl2F2. These two subroutines are called with high frequency and are rewritten with the CUDA interface. Simultaneously, kernel directives and OpenACC directives are adopted to simplify specific operations on the device, and multiple asynchronous streams are also created to overlap function execution and data replication. For the compilation phase, the option -O4 is selected to perform the optimization at the highest level. The option -fast and -fastsse are also utilized to launch the 64-bit Single Instruction Multiple Data (SIMD) instruction and implement cache alignment and flush.In each iteration of PGAPSO, the calculation of the fitness value for each particle is relatively independent. According to its features, multiprocess techniques are suitable for executing these tasks concurrently. Here, MPI is adopted as the parallel framework to accelerate the algorithm. The pseudocode describes the flow of the parallel PGAPSO in Table3.Table 3
The pseudocode of the procedures in PGAPSO.
Algorithm 1 pseudocode for PGAPSO1: transform training data through PCA to obtain principal componentsLd with m dimensions2: initialize population3: FORiter = 1 to iter = itermax DO4:  FOR each processi DO5:   restore solution matrixXi into original space via Xi∗Ld6:   calculate fitness valueJs0=FMt0⟶tS0+P∗Ld under the constraint ϵ7:  END FOR8:  IF norm of particle speed≤ξ THEN9:   select individuals according tops=JsXi′/∑JsX′10:   WHILE count≤ N/2 DO11:    crossover with probabilitypc12:   END WHILE13:   WHILE count≤ N DO14:    mutate with probabilitypm15:   END WHILE16:   compare fitness values between new generation and parent individuals17:  END IF18:  update particle speed viaXpbk and Xgbk19:  update particle position20: END FORThese processes are divided into two groups: one is a master process; the others are slave processes. Figure4 illustrates the parallel architecture of PGAPSO for solving CNOP. At each iteration, the master process allocates calculation tasks to slave processes. For each slave process, perturbation under constraint condition is generated and is superimposed on the basic state of CESM. CESM, which is paralleled with MPI, OpenMP, and CUDA, is called to perform the integration. The fitness values are calculated via equation (9) and are gathered by the master process. After that, the master process broadcasts the optimal global value to slave processes via MPI. If particles’ speed norm is less than the threshold value, the crossover and mutation operations are conducted. Then, the master process compares the current particles and updates the information of optimal particles at the end of each iteration. With the help of MPI, the performance of PGAPSO can be significantly enhanced.Figure 4
The parallel architecture of PGAPSO for solving CNOP.


## 3. Experimental and Results

### 3.1. Experimental Environment

Experiments in this paper are conducted on the Tianhe-2 supercomputer, which is located in the National Supercomputer Center in Guangzhou, China. Each node consists of 2 Intel Ivy Bridge Xeon processors connected by Intel QuickPath Interconnect. NVIDIA Tesla K80 GPUs on Tianhe-2 are used in our GPU-based scheme for CESM acceleration. Each Tesla K80 GPU has 4992 CUDA cores, and its double-precision performance is up to 2.91TFLOPS. Data transmission between CPUs and GPUs depends on PCI-e 3.0 bus with 40 lanes.


### 3.2. Dimensions of Solution Space

The first step is to decompose the principal component from the original sample. The dimension of solution space is selected according to the accumulative variance ratio. Figure5 displays the accumulative variance ratio for different component dimensions.Figure 5
The accumulative variance ratio at different dimensions of principal components.The fitness value and the elapsed steps for different component dimensions are shown in Figure6. Owing to more principal components, the feature space with a larger accumulative variance ratio contains more information of the original space, and a better fitness value can be obtained. However, as the search space becomes larger, it would take a long time or more iterations to find the optimal value. From the figure, it is illustrated that as the dimension of principal components increases, the fitness value increases slowly, and the step number increases obviously. The difference between the growth rate of step and the growth rate of fitness value for dimension = 50 is 9.63%, which is smaller than that of dimension >50, and it has an appropriate fitness value and run time.Figure 6
The fitness value and corresponding steps for different dimensions.Combining these two aspects, the top 50 eigenvectors are selected as the principal components corresponding to the cumulative explained variance ratio of just over 90% to balance the efficiency and effectiveness. In other words, the OPRs are solved in feature space with 50 dimensions. The feature space withm dimensions is far smaller than the original one, and the scale of the original space is 19538 times larger than that of feature space.


### 3.3. The Fusion Threshold

In the PGAPSO, the threshold for GA fusion is determined by the convergence of the PPSO, which is presented by the variation of the fitness value during the iteration process. The search state of the algorithm can be inferred by analyzing the change of the fitness value between two adjacent steps. The algorithm is considered to be at a local optimum when no updating in fitness value for several iterations. A significant change in the fitness value between two adjacent steps indicates that the algorithm jumps out of the local optimum, while the fitness value remains unchanged for a number of steps indicating that the algorithm converges. The variation of the fitness value and the mean speed norm of the single PPSO algorithm with step number is shown in Figure7, where it can be seen that the PPSO tends to fall into a local optimum at around Step 15, and the corresponding mean speed norm is about 20. Thus, the threshold ξ is set to 20 to avoid premature convergence.Figure 7
The fitness value and the mean speed norm within the maximum step.


### 3.4. Simulation Duration for OPR

Since the NAO has an approximated e-folding period of two weeks, the most appropriate simulation duration needs to be determined within 15 days. 5 days, 7 days, 10 days, and 15 days are selected as the optimization times to observe the largest variation of the NAOI amplitude. Using the PGAPSO with 50 iterations, the distribution of the optimal fitness value at different simulation durations is illustrated in Figure8. Each box extends from the lower (Q1) to the upper quartile (Q3), and the line in the box refers to the median. Assume that interquartile range (IQR) = Q3–Q1, then the lower whiskers extend from Q1 to Q1–1.5IQR, and the upper whiskers extend from Q3 to Q3 + 1.5IQR. From the box plot, it can be seen that the 5-day simulation has the narrowest range of the ΔNAOI, especially for NAO−. The 7-day simulation has a wider variation range of values for NAO−, and the 10-day optimization is more inclined to evolve into the NAO−. By contrast, only the 15-day simulation can steadily trigger the NAO events for both two types of phases. The maximum and the minimum values of ΔNAOI are 2.23 and −3.91, which is much larger than 1 or less than -1.Figure 8
The distribution of optimal fitness values at different simulation durations.Then whether the OPRs obtained by PGAPSO can trigger the NAO events is need to be verified. The SLP difference patterns are acquired by subtracting the SLP data of the reference state from the output, and these patterns for different simulation durations, which are generated by 50 iterations, are shown in Figure9. The patterns of 5-day and 7-day simulations are in a state of instability with multiple dispersed pressure centers. It is found that a strong negative center appears in the mode of NAO− in 7-day simulation, and this is the reason for the relatively lower ΔNAOI− in Figure 8. In the same way, the strong positive center in the pattern of 10-day simulation also causes the abnormal ΔNAOI−. The typical feature of the NAO event is the dipole SLP mode located near Iceland and Azores. As shown in Figure 3, in a typical pattern of NAO events, the positive and negative SLP fields are symmetrical along 50° N in the North Atlantic region. The negative SLP field is located in the north when NAO− event occurs, while the positive SLP field is in the north during NAO+ event. Although the approximated structure of NAO events appears in the 7-day and 10-day optimization, the pressure cores are dispersed, and they have not formed a large area of the SLP field with the same phase. As for the 15-day integration, the dipole centers form and migrate across the Atlantic Ocean, which is particularly evident in NAO−. Therefore, the period of 15 days is selected as the simulation duration in our experiments.Figure 9
The difference of SLP patterns between perturbation state and reference state at different simulation durations.


### 3.5. Evolution of NAOI and SLP Pattern

To visualize the effects of PGAPSO, it is compared with several adjoint-free methods, including the breeding of growth mode (BGM), random method, single intelligent algorithms like PPSO, PCAGA, and the PCA-based bacterial foraging algorithm (PBFOA) [60]. The BGM is one of the representative methods of the non-CNOP method, which has been put into operation by the National Centers for Environmental Prediction (NCEP) [61]. Initially, a small perturbation is introduced in the control analysis. After model integration, the difference between the control and perturbation forecasts is reshaped into the size of the initial perturbation, and this difference field is superimposed on the new analysis field. After several iterations, the fastest-growing perturbation is generated. Such iterative process is designed to simulate the generation of growing errors through successive short-term forecasts, as shown in Figure 10. This method selects the fastest-growing perturbation over the forecast period and is used as the baseline method for OPR identification.Figure 10
The forecast procedure of the BGM method.Another procedure is the random method, which is a method based on the CNOP approach without the thought of intelligence. In each iteration of the random method, particles always search with the random position, and the global best position is recorded. The search directions of particles are not affected by the previous iterations and are completely random. In addition, the PCA-based single swarm intelligence algorithms PPSO, PCAGA, and PBFOA are also compared with it. Thereinto, the PBFOA is a global optimization algorithm that solves non-gradient optimization by simulating bacterial chemotaxis in the solution space.Owing to the significance of forecast skills in winter [62], we select two cases (Case 1 and Case 2) in the winter of 53rd model year to simulate the NAO within 15 days. The reference states of these two cases are displayed in Figure 11(a). With the similar initial NAOI but different initial conditions, Case 1 and Case 2 have opposite phases, NAO+ for Case 1 and NAO− for Case 2. Although the NAO event did not occur on Day 15 in these two cases, the final states with the relative higher (lower) NAOI display the trend of evolving into NAO events. This ensures that NAO events in perturbation states are triggered by OPRs. Therefore, Case 1 and Case 2 can be viewed as representative cases to explore the role of OPRs.Figure 11
The trends of the NAOI amplitude for different methods and reference states. (a) Reference state. (b) Case 1. (c) Case 2.
(a)(b)(c)The NAOI amplitude flows for these methods are displayed in Figure11. The black dashed line refers to the reference state, which is the final state acquired without perturbations. Compared with the BGM method, CNOP-type perturbations generally perform better on both NAO+ and NAO−. The perturbations generated by the BGM are effective for forming the structures that reflect the NAO events in NAO−. However, the BGM does not work well in NAO+ whether in Case 1 or Case 2. Besides, the NAOI curves corresponding to the BGM method have a similar varying tendency in the first 10 days and evolve into the contrary phases in the final days of the simulation duration. Since the search directions of random vectors are entirely random, they do not cause the apparent variation. As for the single intelligence algorithms including PCAGA, PPSO, and PBFOA, their results are better than that of the BGM method in NAO+, but approximate to or slightly worse than the BGM method in NAO−. Notably, only PGAPSO and PPSO trigger the extremely high NAOI in Case 2, and PGAPSO outperforms other algorithms in both Case 1 and Case 2. The development of perturbations acquired by the BGM method or random method depends on the basic state to some extent, while the perturbation searched by PGAPSO steadily makes NAOI greater than 2 in both NAO+ and NAO−, and the results are almost not affected by the initial states.Figure12 shows the SLP patterns using these methods in Case 1. For the NAO+, although the center intensity is slightly weak in the pattern of PGAPSO, it presents the principal characteristics of NAO+ events. The opposite pressure fields generate in the center of the Atlantic sector and are symmetrical about the 50°N. In BGM’s pattern, a strong positive core appears in the center of the Atlantic Ocean, while the negative pressure field migrates to the Norwegian Sea, resulting in the overall structure tilts to the east. The pattern of the random method is almost occupied by positive pressure, which is mainly located in Irish islands. The pattern of the PPSO is close to that of the PGAPSO, with lower intensity, while patterns of PCAGA and PBFOA have multiple negative pressure cores and are not in accord with the typical pattern. Unlike the NAO+, the NAO− displays similar structures in these patterns, which have the typical NAO mode with a strong positive center located in Iceland and the negative center(s) around the Azores. The difference is that the patterns of the BGM method, the random method, the PCAGA, and the PBFOA include several discrete negative cores. In Case 1, the OPRs solved by PGAPSO can trigger two types of NAO events with typical dipole mode.Figure 12
The SLP anomaly patterns obtained by multiple methods forNAO+ and NAO− in Case 1.The situation in Case 2 is illustrated in Figure13. For NAO+, the OPR obtained by PGAPSO forms an SLP pattern with noticeable features of the NAO+ event, while the BGM method, the random method, and the PBFOA generate irregular SLP fields with scattered centers. The PPSO has an east-tilt negative pressure center with symmetry, and the pattern of the PCAGA can also be viewed as the NAO event. For NAO−, the intensities are much lower than the NAO+, especially evident in patterns of the PCAGA, the PPSO, and the PBFOA. From Figure 11, it is found that the perturbation state has a larger probability of moving to the contrary phase of the reference. while OPRs obtained by the PGAPSO can still make the basic state evolve into the NAO event that has the same phase as the reference state.Figure 13
The SLP anomaly patterns obtained by multiple methods forNAO+ and NAO− in Case 2.To observe the evolution of the NAO events during the 15-day simulation period, SLP difference fields triggered byCNOPPO and CNOPNE on Day 1, Day 5, Day 10, and Day 15 are plotted in Figure 14. It can be seen from the figure that the SLP difference patterns do not show any changes on the first day—just after the perturbations are superimposed. On Day 5, the difference is still unnoticeable, and a relatively concentrated pressure center appears in the Greenland Sea only in the mode for NAO−. As of Day 10, the change in pressure difference spreads outward from the center of the Arctic, and the area of the negative (positive) pressure field in the Greenland Sea widens, with an increase in intensity. At the prediction time (Day 15), as the NAO+ event occurs, the positive pressure over the Arctic Ocean reveals a regime shift from Day 10, and the negative pressure field around Iceland evolves into a strong one. The positive pressure field on the North Pacific Ocean also has an increasing trend. For the NAO− event, the positive cores around the North Pacific Ocean gather toward the North Pole, and the positive pressure center on Iceland and Greenland gradually expands and reaches its peak. In this process, the most obvious stage of change is from Day 10 to Day 15, and the large difference of pressure and the dipole mode is formed at this stage. This also confirms the conclusion that the role of the nonlinear process is mainly at the end of the prediction phase.Figure 14
The nonlinear evolution forCNOPPO and CNOPNE on Day 1, Day 5, Day 10, and Day 15 at SLP field (Pa) with a simulation time of 15 days.


### 3.6. CNOP Pattern

In this section, the structures of the CNOPs that trigger the NAO event are explored. Figure15 shows the CNOPs obtained using the PGAPSO method for Case 1 in the previous section, including zonal wind U, meridional wind V, and temperature T in the near-surface layer. Figure 15(b) is the CNOP corresponding to the NAO+, and the positive temperature structure appears in the Arctic Ocean region, accompanied by the wind from the Arctic to the Pacific Ocean. A negative temperature field distributes near Iceland, with the wind perturbation direction from the Arctic to the Atlantic. The wind direction of negative-phase CNOP is opposite to that of positive-phase CNOP in the Arctic Ocean and near Iceland, and the intensity of temperature is stronger than that of positive-phase CNOP.Figure 15
The CNOPs forNAO+ (a) and NAO− (b) in Case 1, consisting of the U, V, and T in the 25th layer.
(a)(b)For Case 2, the CNOP of theNAO+ has a higher temperature intensity, while the CNOP of the NAO− has smaller values in temperatures, along with the monotonous distribution characteristics. As can be seen from Figure 11 and Figure 12, under the circumstance of the relative high reference NAOI in Case 1, the amplitude of the CNOP for NAO+ is tiny compared to the NAO−, and the initial state after perturbations superposition is more liable to develop into a NAO− state. The situation in Case 2 is similar to Case 1. The reference NAOI is in a negative phase, and the perturbation that can trigger a strong NAO+ event has a larger intensity. This could mean that the simulation result is more inclined to migrate to the opposite phase as the reference state. It can also be seen from Figures 15 and 16 that, for the two cases whose reference states are widely different, the patterns and value ranges of the CNOPs are also substantially different.Figure 16
The CNOPs forNAO+ (a) and NAO− (b) in Case 2, consisting of the U, V, and T in the 25th layer.
(a)(b)The distribution of the components of specific humidityQ, surface pressure Π, and surface geopotential Φ in the OPRs for Case 1 and Case 2 is shown in Figure 17. The specific humidity has a similar spatial structure height in NAO+ and NAO− for both cases, but with a slightly higher intensity in Case 2. In addition, the NAO+ and NAO− of Case 1 show opposite distributions of specific humidity in the region between the eastern side of Iceland and Scandinavia. The spatial structure of surface pressure differs slightly between NAO+ in Case 2 and the OPRs of other three groups, specifically the Norwegian Sea and the northern North American region centered on the Parry Islands, with only NAO+ of Case 2 being negative in this region, while the other three groups have positive values. For surface geopotential, the OPR component of Case 2 is stronger than that of Case 1, with the common denominator being the positive radial component from the North Pole to the Chulechi Sea. Similar to the surface pressure, the surface geopotential has a similar distribution for both types of NAO events in Case 1 and NAO− in Case 2, with positive values over Greenland, positive values over Iceland, positive values over the Barents Sea coast, etc. In contrast, the surface geopotential of the OPR in Case 2’s NAO+ shows a large-scale negative perturbation.Figure 17
OPRs corresponding toNAO+ and NAO− in Case 1 and Case 2, containing Q, Π, and Φ.Large uncertainties have been triggered by the OPRs containing several variables, and the OPR only with temperature is also found to have a significant impact on the NAO. Following the above procedure, temperature perturbations are limited under a constrained condition of1/D∗∫D∫01T′2dσdD≤100 and are superimposed on the 25th layer, corresponding to the level of 7.389 (near the surface) in the same region with the initial field of Case 1. Using the PGAPSO, the CNOPs composed of temperature are obtained, as illustrated in Figure 18. The patterns of CNOP of NAO+ and CNOP of NAO− have almost converse structures in the North Atlantic sector. There exists a noticeable pressure difference between Greenland and Iceland, with several centers in the mid-to-high latitudes and small cores around the Arctic region. Besides, the positive anomaly in eastern Europe is also conducive to the formation of the dipole. It aligns with the hypotheses that atmospherical temperature gradients will result in anomalous poleward atmospherical heat transport and an increased probability of the NAO reaching its high-index state.Figure 18
The temperature CNOPs superimposing on the Arctic region forNAO+ (a) and NAO− (b).
(a)(b)


### 3.7. Performance Analysis

Figure19 compares the optimal fitness value of NAO+ and NAO− obtained by these methods. Thereinto, the BGM method has good performance in the OPR identification of NAO−, but suffers from deficiencies in the NAO+, while the random method is less effective. Among the single intelligence algorithms, the recognition result of OPR for NAO+ is generally better than that of NAO−’s OPR, with PPSO achieving relatively higher value. The effectiveness of PGAPSO is further improved, especially in the case of NAO−. This also demonstrates the ability of the hybrid algorithm to improve the search for superiority over the single intelligence algorithms.Figure 19
The NAOI values ofNAO+ and NAO− using several methods.Figure20 shows the stability of these algorithms, and the stability is measured using the standard deviation. As can be seen from Figure 10, the results in each step of the BGM are determined by the previous step; thus the result is a constant with a standard deviation of 0. As for other methods, the random method is the least stable; PPSO is more stable for OPR identification at NAO− but fluctuates greatly at NAO+; PCAGA, PBFOA, and PGAPSO have similar stability, with PGAPSO having the highest stability at NAO+. It is proved that the robustness of PGAPSO is also reasonable.Figure 20
The stability of different methods (the standard deviation of the BGM method is 0).In this paper, the convergence of the above algorithm is explored by recording the changes of the fitness values, as shown in Figure21. The convergence of the BGM method is not analyzed since its results are constants. As mentioned in Figure 7, if the fitness values continue to remain unchanged for a long time, the algorithm is considered to be converged. The random method converges faster but with lower fitness values due to the lack of guidance on the search direction, resulting in lower efficiency; PCAGA continues to be at a local optimum around 19 steps, resulting in a limited effect of the algorithm; PPSO jumps out of the local optimum more often and has a stronger global search capability; PGAPSO still achieves better fitness values in NAO+ and NAO−. In the early stage of the algorithm’s search, PGAPSO has the most times to jump out of the local optimum, which led to frequent updates of the optimal fitness values and achieves better results, converging at 35 steps (NAO+) and 32 steps (NAO−), respectively. PGAPSO has an increased update number of optimal fitness values compared to PPSO, indicating that the randomness of the algorithm is improved by combining the GA. The capacity for jumping out of the local optimum is enhanced and premature convergence is avoided. The convergence of the PBFOA performs between PCAGA and PPSO. The results show that the PGAPSO algorithm has better global search capability in early stage and local search capability in late stage, resulting in strong convergence (Figure 21).Figure 21
The variation of NAOI for the abovementioned methods in 50 steps.To manifest the performance improvement of parallel PGAPSO, Figure22 compares the runtime of parallel PGAPSO and serial PGAPSO in one iteration. Due to the limitation of computing resources, experiments are conducted on no more than 1080 CPU cores. From Figure 22, it can be seen that when the number of CPU cores is more than 720, it will take longer to run the serial algorithm. When the overmuch CPU cores are assigned to the serial program, the frequent communications would make the runtime increase. In contrast, the parallel version can make full use of resources, and its execution time keeps going down.Figure 22
The execution time of serial PGAPSO and parallel PGAPSO with different number of CPU cores.The speedup ratio, which raises with the increasing CPU cores’ number, is displayed in Table4. With 1080 CPU cores, PGAPSO, based on the parallel scheme, achieves a speedup of 40.2 × compared to its serial version.Table 4
The speedup of parallel PGAPSO compared with serial PGAPSO.
Number of CPU coresSpeedup ratio24016.548017.972018.784026.396033.2108040.2


## 3.1. Experimental Environment

Experiments in this paper are conducted on the Tianhe-2 supercomputer, which is located in the National Supercomputer Center in Guangzhou, China. Each node consists of 2 Intel Ivy Bridge Xeon processors connected by Intel QuickPath Interconnect. NVIDIA Tesla K80 GPUs on Tianhe-2 are used in our GPU-based scheme for CESM acceleration. Each Tesla K80 GPU has 4992 CUDA cores, and its double-precision performance is up to 2.91TFLOPS. Data transmission between CPUs and GPUs depends on PCI-e 3.0 bus with 40 lanes.


## 3.2. Dimensions of Solution Space

The first step is to decompose the principal component from the original sample. The dimension of solution space is selected according to the accumulative variance ratio. Figure5 displays the accumulative variance ratio for different component dimensions.Figure 5
The accumulative variance ratio at different dimensions of principal components.The fitness value and the elapsed steps for different component dimensions are shown in Figure6. Owing to more principal components, the feature space with a larger accumulative variance ratio contains more information of the original space, and a better fitness value can be obtained. However, as the search space becomes larger, it would take a long time or more iterations to find the optimal value. From the figure, it is illustrated that as the dimension of principal components increases, the fitness value increases slowly, and the step number increases obviously. The difference between the growth rate of step and the growth rate of fitness value for dimension = 50 is 9.63%, which is smaller than that of dimension >50, and it has an appropriate fitness value and run time.Figure 6
The fitness value and corresponding steps for different dimensions.Combining these two aspects, the top 50 eigenvectors are selected as the principal components corresponding to the cumulative explained variance ratio of just over 90% to balance the efficiency and effectiveness. In other words, the OPRs are solved in feature space with 50 dimensions. The feature space withm dimensions is far smaller than the original one, and the scale of the original space is 19538 times larger than that of feature space.


## 3.3. The Fusion Threshold

In the PGAPSO, the threshold for GA fusion is determined by the convergence of the PPSO, which is presented by the variation of the fitness value during the iteration process. The search state of the algorithm can be inferred by analyzing the change of the fitness value between two adjacent steps. The algorithm is considered to be at a local optimum when no updating in fitness value for several iterations. A significant change in the fitness value between two adjacent steps indicates that the algorithm jumps out of the local optimum, while the fitness value remains unchanged for a number of steps indicating that the algorithm converges. The variation of the fitness value and the mean speed norm of the single PPSO algorithm with step number is shown in Figure7, where it can be seen that the PPSO tends to fall into a local optimum at around Step 15, and the corresponding mean speed norm is about 20. Thus, the threshold ξ is set to 20 to avoid premature convergence.Figure 7
The fitness value and the mean speed norm within the maximum step.


## 3.4. Simulation Duration for OPR

Since the NAO has an approximated e-folding period of two weeks, the most appropriate simulation duration needs to be determined within 15 days. 5 days, 7 days, 10 days, and 15 days are selected as the optimization times to observe the largest variation of the NAOI amplitude. Using the PGAPSO with 50 iterations, the distribution of the optimal fitness value at different simulation durations is illustrated in Figure8. Each box extends from the lower (Q1) to the upper quartile (Q3), and the line in the box refers to the median. Assume that interquartile range (IQR) = Q3–Q1, then the lower whiskers extend from Q1 to Q1–1.5IQR, and the upper whiskers extend from Q3 to Q3 + 1.5IQR. From the box plot, it can be seen that the 5-day simulation has the narrowest range of the ΔNAOI, especially for NAO−. The 7-day simulation has a wider variation range of values for NAO−, and the 10-day optimization is more inclined to evolve into the NAO−. By contrast, only the 15-day simulation can steadily trigger the NAO events for both two types of phases. The maximum and the minimum values of ΔNAOI are 2.23 and −3.91, which is much larger than 1 or less than -1.Figure 8
The distribution of optimal fitness values at different simulation durations.Then whether the OPRs obtained by PGAPSO can trigger the NAO events is need to be verified. The SLP difference patterns are acquired by subtracting the SLP data of the reference state from the output, and these patterns for different simulation durations, which are generated by 50 iterations, are shown in Figure9. The patterns of 5-day and 7-day simulations are in a state of instability with multiple dispersed pressure centers. It is found that a strong negative center appears in the mode of NAO− in 7-day simulation, and this is the reason for the relatively lower ΔNAOI− in Figure 8. In the same way, the strong positive center in the pattern of 10-day simulation also causes the abnormal ΔNAOI−. The typical feature of the NAO event is the dipole SLP mode located near Iceland and Azores. As shown in Figure 3, in a typical pattern of NAO events, the positive and negative SLP fields are symmetrical along 50° N in the North Atlantic region. The negative SLP field is located in the north when NAO− event occurs, while the positive SLP field is in the north during NAO+ event. Although the approximated structure of NAO events appears in the 7-day and 10-day optimization, the pressure cores are dispersed, and they have not formed a large area of the SLP field with the same phase. As for the 15-day integration, the dipole centers form and migrate across the Atlantic Ocean, which is particularly evident in NAO−. Therefore, the period of 15 days is selected as the simulation duration in our experiments.Figure 9
The difference of SLP patterns between perturbation state and reference state at different simulation durations.


## 3.5. Evolution of NAOI and SLP Pattern

To visualize the effects of PGAPSO, it is compared with several adjoint-free methods, including the breeding of growth mode (BGM), random method, single intelligent algorithms like PPSO, PCAGA, and the PCA-based bacterial foraging algorithm (PBFOA) [60]. The BGM is one of the representative methods of the non-CNOP method, which has been put into operation by the National Centers for Environmental Prediction (NCEP) [61]. Initially, a small perturbation is introduced in the control analysis. After model integration, the difference between the control and perturbation forecasts is reshaped into the size of the initial perturbation, and this difference field is superimposed on the new analysis field. After several iterations, the fastest-growing perturbation is generated. Such iterative process is designed to simulate the generation of growing errors through successive short-term forecasts, as shown in Figure 10. This method selects the fastest-growing perturbation over the forecast period and is used as the baseline method for OPR identification.Figure 10
The forecast procedure of the BGM method.Another procedure is the random method, which is a method based on the CNOP approach without the thought of intelligence. In each iteration of the random method, particles always search with the random position, and the global best position is recorded. The search directions of particles are not affected by the previous iterations and are completely random. In addition, the PCA-based single swarm intelligence algorithms PPSO, PCAGA, and PBFOA are also compared with it. Thereinto, the PBFOA is a global optimization algorithm that solves non-gradient optimization by simulating bacterial chemotaxis in the solution space.Owing to the significance of forecast skills in winter [62], we select two cases (Case 1 and Case 2) in the winter of 53rd model year to simulate the NAO within 15 days. The reference states of these two cases are displayed in Figure 11(a). With the similar initial NAOI but different initial conditions, Case 1 and Case 2 have opposite phases, NAO+ for Case 1 and NAO− for Case 2. Although the NAO event did not occur on Day 15 in these two cases, the final states with the relative higher (lower) NAOI display the trend of evolving into NAO events. This ensures that NAO events in perturbation states are triggered by OPRs. Therefore, Case 1 and Case 2 can be viewed as representative cases to explore the role of OPRs.Figure 11
The trends of the NAOI amplitude for different methods and reference states. (a) Reference state. (b) Case 1. (c) Case 2.
(a)(b)(c)The NAOI amplitude flows for these methods are displayed in Figure11. The black dashed line refers to the reference state, which is the final state acquired without perturbations. Compared with the BGM method, CNOP-type perturbations generally perform better on both NAO+ and NAO−. The perturbations generated by the BGM are effective for forming the structures that reflect the NAO events in NAO−. However, the BGM does not work well in NAO+ whether in Case 1 or Case 2. Besides, the NAOI curves corresponding to the BGM method have a similar varying tendency in the first 10 days and evolve into the contrary phases in the final days of the simulation duration. Since the search directions of random vectors are entirely random, they do not cause the apparent variation. As for the single intelligence algorithms including PCAGA, PPSO, and PBFOA, their results are better than that of the BGM method in NAO+, but approximate to or slightly worse than the BGM method in NAO−. Notably, only PGAPSO and PPSO trigger the extremely high NAOI in Case 2, and PGAPSO outperforms other algorithms in both Case 1 and Case 2. The development of perturbations acquired by the BGM method or random method depends on the basic state to some extent, while the perturbation searched by PGAPSO steadily makes NAOI greater than 2 in both NAO+ and NAO−, and the results are almost not affected by the initial states.Figure12 shows the SLP patterns using these methods in Case 1. For the NAO+, although the center intensity is slightly weak in the pattern of PGAPSO, it presents the principal characteristics of NAO+ events. The opposite pressure fields generate in the center of the Atlantic sector and are symmetrical about the 50°N. In BGM’s pattern, a strong positive core appears in the center of the Atlantic Ocean, while the negative pressure field migrates to the Norwegian Sea, resulting in the overall structure tilts to the east. The pattern of the random method is almost occupied by positive pressure, which is mainly located in Irish islands. The pattern of the PPSO is close to that of the PGAPSO, with lower intensity, while patterns of PCAGA and PBFOA have multiple negative pressure cores and are not in accord with the typical pattern. Unlike the NAO+, the NAO− displays similar structures in these patterns, which have the typical NAO mode with a strong positive center located in Iceland and the negative center(s) around the Azores. The difference is that the patterns of the BGM method, the random method, the PCAGA, and the PBFOA include several discrete negative cores. In Case 1, the OPRs solved by PGAPSO can trigger two types of NAO events with typical dipole mode.Figure 12
The SLP anomaly patterns obtained by multiple methods forNAO+ and NAO− in Case 1.The situation in Case 2 is illustrated in Figure13. For NAO+, the OPR obtained by PGAPSO forms an SLP pattern with noticeable features of the NAO+ event, while the BGM method, the random method, and the PBFOA generate irregular SLP fields with scattered centers. The PPSO has an east-tilt negative pressure center with symmetry, and the pattern of the PCAGA can also be viewed as the NAO event. For NAO−, the intensities are much lower than the NAO+, especially evident in patterns of the PCAGA, the PPSO, and the PBFOA. From Figure 11, it is found that the perturbation state has a larger probability of moving to the contrary phase of the reference. while OPRs obtained by the PGAPSO can still make the basic state evolve into the NAO event that has the same phase as the reference state.Figure 13
The SLP anomaly patterns obtained by multiple methods forNAO+ and NAO− in Case 2.To observe the evolution of the NAO events during the 15-day simulation period, SLP difference fields triggered byCNOPPO and CNOPNE on Day 1, Day 5, Day 10, and Day 15 are plotted in Figure 14. It can be seen from the figure that the SLP difference patterns do not show any changes on the first day—just after the perturbations are superimposed. On Day 5, the difference is still unnoticeable, and a relatively concentrated pressure center appears in the Greenland Sea only in the mode for NAO−. As of Day 10, the change in pressure difference spreads outward from the center of the Arctic, and the area of the negative (positive) pressure field in the Greenland Sea widens, with an increase in intensity. At the prediction time (Day 15), as the NAO+ event occurs, the positive pressure over the Arctic Ocean reveals a regime shift from Day 10, and the negative pressure field around Iceland evolves into a strong one. The positive pressure field on the North Pacific Ocean also has an increasing trend. For the NAO− event, the positive cores around the North Pacific Ocean gather toward the North Pole, and the positive pressure center on Iceland and Greenland gradually expands and reaches its peak. In this process, the most obvious stage of change is from Day 10 to Day 15, and the large difference of pressure and the dipole mode is formed at this stage. This also confirms the conclusion that the role of the nonlinear process is mainly at the end of the prediction phase.Figure 14
The nonlinear evolution forCNOPPO and CNOPNE on Day 1, Day 5, Day 10, and Day 15 at SLP field (Pa) with a simulation time of 15 days.


## 3.6. CNOP Pattern

In this section, the structures of the CNOPs that trigger the NAO event are explored. Figure15 shows the CNOPs obtained using the PGAPSO method for Case 1 in the previous section, including zonal wind U, meridional wind V, and temperature T in the near-surface layer. Figure 15(b) is the CNOP corresponding to the NAO+, and the positive temperature structure appears in the Arctic Ocean region, accompanied by the wind from the Arctic to the Pacific Ocean. A negative temperature field distributes near Iceland, with the wind perturbation direction from the Arctic to the Atlantic. The wind direction of negative-phase CNOP is opposite to that of positive-phase CNOP in the Arctic Ocean and near Iceland, and the intensity of temperature is stronger than that of positive-phase CNOP.Figure 15
The CNOPs forNAO+ (a) and NAO− (b) in Case 1, consisting of the U, V, and T in the 25th layer.
(a)(b)For Case 2, the CNOP of theNAO+ has a higher temperature intensity, while the CNOP of the NAO− has smaller values in temperatures, along with the monotonous distribution characteristics. As can be seen from Figure 11 and Figure 12, under the circumstance of the relative high reference NAOI in Case 1, the amplitude of the CNOP for NAO+ is tiny compared to the NAO−, and the initial state after perturbations superposition is more liable to develop into a NAO− state. The situation in Case 2 is similar to Case 1. The reference NAOI is in a negative phase, and the perturbation that can trigger a strong NAO+ event has a larger intensity. This could mean that the simulation result is more inclined to migrate to the opposite phase as the reference state. It can also be seen from Figures 15 and 16 that, for the two cases whose reference states are widely different, the patterns and value ranges of the CNOPs are also substantially different.Figure 16
The CNOPs forNAO+ (a) and NAO− (b) in Case 2, consisting of the U, V, and T in the 25th layer.
(a)(b)The distribution of the components of specific humidityQ, surface pressure Π, and surface geopotential Φ in the OPRs for Case 1 and Case 2 is shown in Figure 17. The specific humidity has a similar spatial structure height in NAO+ and NAO− for both cases, but with a slightly higher intensity in Case 2. In addition, the NAO+ and NAO− of Case 1 show opposite distributions of specific humidity in the region between the eastern side of Iceland and Scandinavia. The spatial structure of surface pressure differs slightly between NAO+ in Case 2 and the OPRs of other three groups, specifically the Norwegian Sea and the northern North American region centered on the Parry Islands, with only NAO+ of Case 2 being negative in this region, while the other three groups have positive values. For surface geopotential, the OPR component of Case 2 is stronger than that of Case 1, with the common denominator being the positive radial component from the North Pole to the Chulechi Sea. Similar to the surface pressure, the surface geopotential has a similar distribution for both types of NAO events in Case 1 and NAO− in Case 2, with positive values over Greenland, positive values over Iceland, positive values over the Barents Sea coast, etc. In contrast, the surface geopotential of the OPR in Case 2’s NAO+ shows a large-scale negative perturbation.Figure 17
OPRs corresponding toNAO+ and NAO− in Case 1 and Case 2, containing Q, Π, and Φ.Large uncertainties have been triggered by the OPRs containing several variables, and the OPR only with temperature is also found to have a significant impact on the NAO. Following the above procedure, temperature perturbations are limited under a constrained condition of1/D∗∫D∫01T′2dσdD≤100 and are superimposed on the 25th layer, corresponding to the level of 7.389 (near the surface) in the same region with the initial field of Case 1. Using the PGAPSO, the CNOPs composed of temperature are obtained, as illustrated in Figure 18. The patterns of CNOP of NAO+ and CNOP of NAO− have almost converse structures in the North Atlantic sector. There exists a noticeable pressure difference between Greenland and Iceland, with several centers in the mid-to-high latitudes and small cores around the Arctic region. Besides, the positive anomaly in eastern Europe is also conducive to the formation of the dipole. It aligns with the hypotheses that atmospherical temperature gradients will result in anomalous poleward atmospherical heat transport and an increased probability of the NAO reaching its high-index state.Figure 18
The temperature CNOPs superimposing on the Arctic region forNAO+ (a) and NAO− (b).
(a)(b)


## 3.7. Performance Analysis

Figure19 compares the optimal fitness value of NAO+ and NAO− obtained by these methods. Thereinto, the BGM method has good performance in the OPR identification of NAO−, but suffers from deficiencies in the NAO+, while the random method is less effective. Among the single intelligence algorithms, the recognition result of OPR for NAO+ is generally better than that of NAO−’s OPR, with PPSO achieving relatively higher value. The effectiveness of PGAPSO is further improved, especially in the case of NAO−. This also demonstrates the ability of the hybrid algorithm to improve the search for superiority over the single intelligence algorithms.Figure 19
The NAOI values ofNAO+ and NAO− using several methods.Figure20 shows the stability of these algorithms, and the stability is measured using the standard deviation. As can be seen from Figure 10, the results in each step of the BGM are determined by the previous step; thus the result is a constant with a standard deviation of 0. As for other methods, the random method is the least stable; PPSO is more stable for OPR identification at NAO− but fluctuates greatly at NAO+; PCAGA, PBFOA, and PGAPSO have similar stability, with PGAPSO having the highest stability at NAO+. It is proved that the robustness of PGAPSO is also reasonable.Figure 20
The stability of different methods (the standard deviation of the BGM method is 0).In this paper, the convergence of the above algorithm is explored by recording the changes of the fitness values, as shown in Figure21. The convergence of the BGM method is not analyzed since its results are constants. As mentioned in Figure 7, if the fitness values continue to remain unchanged for a long time, the algorithm is considered to be converged. The random method converges faster but with lower fitness values due to the lack of guidance on the search direction, resulting in lower efficiency; PCAGA continues to be at a local optimum around 19 steps, resulting in a limited effect of the algorithm; PPSO jumps out of the local optimum more often and has a stronger global search capability; PGAPSO still achieves better fitness values in NAO+ and NAO−. In the early stage of the algorithm’s search, PGAPSO has the most times to jump out of the local optimum, which led to frequent updates of the optimal fitness values and achieves better results, converging at 35 steps (NAO+) and 32 steps (NAO−), respectively. PGAPSO has an increased update number of optimal fitness values compared to PPSO, indicating that the randomness of the algorithm is improved by combining the GA. The capacity for jumping out of the local optimum is enhanced and premature convergence is avoided. The convergence of the PBFOA performs between PCAGA and PPSO. The results show that the PGAPSO algorithm has better global search capability in early stage and local search capability in late stage, resulting in strong convergence (Figure 21).Figure 21
The variation of NAOI for the abovementioned methods in 50 steps.To manifest the performance improvement of parallel PGAPSO, Figure22 compares the runtime of parallel PGAPSO and serial PGAPSO in one iteration. Due to the limitation of computing resources, experiments are conducted on no more than 1080 CPU cores. From Figure 22, it can be seen that when the number of CPU cores is more than 720, it will take longer to run the serial algorithm. When the overmuch CPU cores are assigned to the serial program, the frequent communications would make the runtime increase. In contrast, the parallel version can make full use of resources, and its execution time keeps going down.Figure 22
The execution time of serial PGAPSO and parallel PGAPSO with different number of CPU cores.The speedup ratio, which raises with the increasing CPU cores’ number, is displayed in Table4. With 1080 CPU cores, PGAPSO, based on the parallel scheme, achieves a speedup of 40.2 × compared to its serial version.Table 4
The speedup of parallel PGAPSO compared with serial PGAPSO.
Number of CPU coresSpeedup ratio24016.548017.972018.784026.396033.2108040.2


## 4. Discussions and Conclusions

Initial condition errors are critical factors that result in uncertainty when simulating and predicting the NAO, and the NAO simulation can be improved by reducing the errors in the initial condition. As a type of initial condition error, the OPR would cause the largest prediction error and eventually evolve into climate events. Therefore, the research of OPRs would help to reveal the dynamic processes related to the NAO events and improve the prediction accuracy.In this paper, a CNOP-based approach PGAPSO is adopted to study OPRs of the NAO using CESM. Since the CESM does not have a corresponding adjoint model, CNOP cannot be solved through the adjoint-based method mentioned in the previous works, such as SQP and SPG2. First, the optimization time is determined from experiments. Then, the perturbations, which contain the variables of zonal wind, meridional wind, temperature, specific humidity, surface pressure, and surface geopotential, are superimposed on the basic state over the Arctic region. The OPRs of two cases in winter (DJF) are investigated with different initial states. To validate the effectiveness of the PGAPSO, the trends of the NAOI amplitude and SLP anomaly patterns are compared to the BGM method, the random method, and single swarm intelligence algorithms like PCAGA, PPSO, and PBFOA. It is indicated that the OPRs obtained by PGAPSO can trigger the NAO events with the typical dipole pattern and have the largestΔNAOI. The SLP variation in the northern hemisphere is traced during simulation time, and the features of teleconnection patterns are identified via the SLP difference mode.The OPRs' structures for bothNAO+ and NAO− in these two cases are also analyzed. The wind directions of NAO+ and NAO− present the opposite mode around the Arctic Ocean and Iceland. The temperature perturbations over Greenland island would promote the occurrence of NAO events; the experiments with OPRs containing the only temperature also confirm this. Specific humidity and surface pressure present a similar spatial structure in the two types of NAO in both cases, while surface potential heights are distinctive in the region from the North Pole to the Chukchi Sea. It is also demonstrated that even slight errors (see NAO− in Case 2) may cause a large uncertainty in simulation.The improvement in solution quality, robustness, and solution efficiency of the PGAPSO algorithm is analyzed in terms of effectiveness, stability, and convergence. Compared with the BGM method, the random method, and single intelligence methods mentioned above, it is confirmed that the improvement in the PGAPSO algorithm helps to improve the quality of the search and the convergence speed and also has a better stability. Moreover, multiple parallel frameworks are applied in this work to improve efficiency. The parallelization mainly consists of two parts: parallelization of the algorithm with MPI and acceleration of CESM using CUDA. It significantly enhances the performance and achieves a speedup of 40.2×.Our future work is to apply the PGAPSO algorithm to study other climatological phenomena with the CNOP method. In this work, the CESM is regarded as a black box program. It is convenient to transplant the solver framework to other numerical models. Our approach will be applied to models that have high dimensions or have no corresponding adjoint model.



---

*Source: 3215039-2022-02-08.xml*