Climate change is one of the key topics in climate science. However, previous research has predominantly concentrated on changes in mean values, and few research examines changes in Probability Distribution Function (PDF). In this study, a novel method called Wasserstein Stability Analysis (WSA) is developed to identify PDF changes, especially the extreme event shift and non-linear physical value constraint variation in climate change. WSA is applied to 21st-century warming slowdown period and is compared with traditional mean-value trend analysis. The result indicates that despite no significant trend, the central-eastern Pacific experienced a decline in hot extremes and an increase in cold extremes, indicating a La Nina-like temperature shift. Further analysis at two Arctic locations suggests sea ice severely restricts the hot extremes of surface air temperature. This impact is diminishing as sea ice melts. Overall, based on detecting PDF changes, WSA is a useful method for re-discovering climate change. How climate system evolves with time is one of the essential topics in climate science. The Intergovernmental Panel on (IPCC AR6) pointed out that the global mean of surface air temperature (SAT) follows a general warming trend since 20th century, which is mainly attributed to the human influence. Consequently, considerable work is devoted to discussing the mean temperature warming trend and its spatial distribution (Delworth & Knutson, 2000;Delworth et al., 2015;Huang et al., 2017;Li et al., 2015). The mean value is one of the most significant statistical characteristics of the data set, as it represents the overall energy condition and even the low-frequency variability of the climate system. However, the human-induced mean value warming of SAT may be partially countered by internal variability or external forcing, thus the warming trend is not always strong and rising at the same pace. As a result, the mean value change, especially the trend based on mean value, cannot give sufficient information about climate change. In the first decade of the 21st century, global warming decelerates significantly (Modak & Mauritsen, 2021;Lee et al., 2015;Delworth et al., 2015;England et al., 2014;Li et al., 2015). However, during the period, a series of climate change are still recorded, such as extreme events (Johnson et al., 2018) and La-Nina-like Pacific temperature pattern (Kosaka & Xie, 2013). Unfortunately, it is difficult to detect these climate change signals by trend or mean value analysis since the relevant signals during the warming slowdown period are quite weak. Recently, a new mathematical tool from optimal transport, called Wasserstein distance (W-distance), sheds more light on this question. W-distance W p (, ) on R d is induced by the study of how to transport mass from a probability distribution  to another distribution  in the cheapest way, which is given by |x  y| p d(x, y) where p  [1, ) and   (, ) means  is a joint distribution of  and . W-distance satisfies the metric axioms and can quantify the distance and similarities between two probability distributions (Figalli & Glaudo, 2021). More details about W-distance could be found in supporting information materials. In this work, we use earth mover's distance, i.e. p = 1. Wasserstein distance has been applied to climate science, such as model evaluation (Vissio et al., 2020), oceanographic data analysis (Hyun et al., 2022), and data assimilation (Tamang, 2022). However, it is rarely applied in the physical analysis of climate data. Based on W-distance, we develop a novel method, named as Wasserstein Stability Analysis (WSA), to discover PDF variability in climate change, which includes signals such as extremes and physical value constraints. Compared to mean value analysis, WSA could help researchers directly identify the significant PDF variability and gives a better insight into climate change. The paper is organized as follows. In section 2, we give a description to the data set and WSA algorithm. Optimal mass transportation and Wasserstein distance are also briefly introduced. In section 3, we apply WSA to global warming slowdown period. Through the new method, we obtain a clear equatorial eastern Pacific signal that is highlighted in previous studies and an evident physical value constraint change in the Arctic that is not fully explored. Finally, in section 4, we draw a conclusion and discuss the future work. The surface 2 meters air temperature data in this work are obtained from ERA5 dataset (0.25 o x 0.25 o , single layer, half day time interval) from 1998 to 2012 (Hersbach et al., 2020). And the spatial range is global. To better represent the large-scale signal and accelerate the data processing, we regrid the original data set into daily data with a horizontal resolution of 2.5 o x 2.5 o by spatial and temporal average. The sea ice data used in the study come from and Sea Surface Temperature monthly data set (HadISST), with 1 o x 1 o horizontal resolution (Rayner et al., 2003). The time range is from 1998 to 2012 and the spatial range is still global. In the null hypothesis of most statistical analyses, we usually assume that the PDF of one variable is stable and follows the same distribution. Following this logic, we claim that the variable has an unstable PDF if the W-distance between its current state and a reference state is out of a specific confidence interval, representing a significant change from one climate stage to another. Therefore, the significant PDF change signal is detected by testing W-distance, and the climate change details can be further explored. Hence, in this part, we develop a novel method named as Wasserstein Stability Analysis (WSA) to evaluate the magnitude of PDF change quantitatively, i.e., to detect unstable PDF change signals. The new method is mainly divided into two steps: W-distance test and PDF analysis. For more information about W-distance test, interested readers could refer to Algorithm 1 in the supporting information material. The details are as follows: (1). W-distance test. Given two samples X and Y , both series are first normalized individually by the following scalar where  is the mean value and  is the standard deviation.
Then Wasserstein distance W 1 (X, Y ) between X and Y are obtained after normalized series. Note that normalization is necessary since the geographic difference may lead to a large variation in W-distance. For example, the high surface temperature variation at high latitudes usually leads to a larger W-distance than its counterpart in tropical zones. Nevertheless, it does not necessarily mean the PDF change at high latitudes is stronger than that at low latitudes. Our significance test is based on the Monte Carro test (Xie et al., 2017). The null hypothesis H 0 of WSA is that the PDF discrepancy between X and Y can be explained by the white noise. We generate two new samples, X N and Y N , by adding white noise to X and Y . Then one could obtain a new W-distance W 1 (X N , Y N ). Such operation is performed 500 times. Then one could get a confidence interval (C.I.) of W-distances. Under the null hypothesis with confidence level 1 If not, H 0 is rejected with the confidence level of 1  . The  in this study has been set as 0.01 and we perform the W-distance test for all locations for spatial pattern. (2). PDF analysis. After detecting the significant PDF change zones where the Wdistance is significant, we plot two PDFs and the corresponding change for each significant change zone, and then explore the climate change mechanism that makes W-distance significant. The next section will use two examples to demonstrate how it functions. During 1998 -2012 period, global mean surface temperature experienced a slower increase, which is referred to as "global warming hiatus" in early studies (Kosaka & Xie, 2013). Recent research indicated that there was no true "hiatus" during that time, but the global surface temperatures were still rising, albeit more slowly (Simmons et al., 2017). This phenomenon suggests that physical processes other than the anthropogenic warming trend play a role in climate change, and the detection of PDF change may provide us with new information for this period. Therefore, by comparing the performance of the WSA with traditional trend analysis, we can evaluate the similarity and difference between these two approaches, which could advance our understanding of climate science and inform future research in this area. has the greatest W-distance, with values ranging from 0.08 to 0.24; nevertheless, there is no discernible trend . Another noticeable W-distance signal could be seen over the Barents-Kara Sea, which is located inside a significant warming trend zone in the Arctic . The following analysis will mainly focus on signals over central eastern Pacific and Barents-Kara Sea to highlight the information that W-distance can provide about climate change. In practice, we found that due to the presence of sea ice, the surface temperature distribution shows a large discontinuous in the Arctic region. As a result, the areal average in this region may smash some information from individual grids. To avoid this, we selected three special sites to do further PDF analysis: 1. Site A (90 o W, 0): one site with a significant W-distance signal but without a significant linear trend; 2. Sites B (66 o E,79 o N) and C (131 o E,79 o N): both sites exhibit significant linear trends. However, Site B has a significant W-distance while C does not. However, in the second period, the hot extremes (more than two standard deviations, roughly 27 o C) decrease by more than ten days per year, while the mild hot events (near one standard deviation) increase by roughly 30 days per year. Meanwhile, the cold extremes (less than -2 standard deviations, 19 o C) increase by more than ten days per year, while the mild cold events (near -1 standard deviation) decrease by more than 30 days per year. The hot extreme decreasing and cold extreme increasing indicate a La Nina-like SAT shift. Many early studies argue that the 21st slowdown period is in the negative phase of Pacific Decadal Oscillation, and this strong low-frequency variability canceled out some parts of the global warming trend (Douville et al., 2015;Kosaka & Xie, 2013;Lee et al., 2015;England et al., 2014). Our W-distance result agrees well with this conclusion and extends the results from the perspective of PDF change. Indeed, even though the mean value does not show a significant trend, we can still observe the La Nina-like SAT shift from extreme event change, which is reflected by the strong PDF change and the corresponding significant W-distance. Site B is another region with a significant W-distance, suggesting a strong PDF change from the 1998-2005 period to the 2005-2012 period and (c)). At Site B, SAT shows significant warming from -10.2 o C to -6.5 o C. Meanwhile, the standard deviation also decreases, indicating an elimination of SAT variation. Since Site B is located within a zone covered by sea ice, the maximum SAT is strongly constrained to being near the frozen line in both periods, and thus this causes a sharp frequency jump near the frozen line ). However, as the mean value increases and the standard deviation decreases during the second period, the frozen line at Site B shifts significantly to the left. As a result, the frequency peak of hot events shifts to the left. Meanwhile, the slope near the frozen line becomes flatter, leading to a smoother PDF shape. In contrast, Site C exhibits a change that follows a similar pattern but is considerably smaller in extent in terms of mean value, standard deviation, and frozen line shift and (d)). Consequently, W-distance at Site C is smaller than that at Site B and fails the 99% significant test. Actually, the existence of sea ice strongly limits the maximum surface temperature and the corresponding surface air temperature near the frozen line. Therefore, the sharp peak shape of SAT PDF appears in regions covered by sea ice. Nevertheless, regions with normal ocean conditions are expected to have a more Gaussian distribution.
Therefore, the alteration in Site B and Site C suggests that the variability in sea ice concentration may play a role in the SAT PDF change from the first to the second period. By checking the full-time series in two periods, we can better understand the PDF change at Site B and Site C. From 1998-2005 period to 2005-2012 period, the interval with +/-1 standard deviation shrinks and moves upward at Site B because of the change in mean value and standard deviation ). Clearly, this change leads to most of the SAT variabilities being below +1 standard deviation and also more cold extremes. For Site C, we can see a similar change but in a much smaller magnitude ). As discussed, this difference can be attributed to the sea ice and its strong temperature constraint. Figure 4 (c) shows the sea ice concentration change between 1998-2005 and 2005-2012. We can see that Site B is in a zone where sea ice is melting the most. Hence the PDF change at Site B is the most intense. As a comparison, Site C also experiences a sea ice loss but to a relatively small extent. As the temperature increases, the sea ice concentration decreases, and the effect of sea ice restriction gradually fades. This is why Site C also displays a similar PDF change but with a less significant level. In summary, the significant W-distance change at Site B implies the removal of physical temperature constraint imposed by sea ice that causes the sharp peak near the frozen line in PDF. Actually, the sea ice-induced PDF change we discuss here is also related to modelling bias. IPCC AR6 indicated that observational-based estimates of global mean surface temperature face a significant issue in areas where sea ice melts or grows, resulting in a switch between air temperature and sea surface temperature. This bias towards reduced warming in anomalies primarily affects analyses combining surface -29.9 -20.0 -10.2 -0.3 9.5 -21.4 -13.9 -6.5 0.9 8.4 0 air temperature anomalies over land and ice with sea surface temperature anomalies over ocean. According to Richardson et al. (2018), this underestimated warming in historical model simulations amounts to about 3% of the observed warming. Given that WSA is able to detect this sea ice-induced change in PDF, it may be useful for future model development and data correction. Overall, by detecting changes in the probability density function (PDF) between the early and later periods during slowdown period, WSA is able to identify the two physical processes during the warming slowdown period: the La Nina-like tropical Pacific temperature shift and physical value constraint due to sea ice. In this work, we develop a novel method called WSA to detect climate change signals. The results indicate that WSA shows some advantages compared with the classic mean value analysis in certain cases. The classic trend analysis based on mean value cannot detect the extreme event signals in the central eastern Pacific Ocean during warming hiatus, while WSA can easily identify the La Nina-like temperature PDF shift. Additionally, due to the strong temperature constraint of sea ice, the sea ice cover area shows a different PDF pattern compared to the open ocean. As a result, unlike mean value analysis detecting significant trend regions, WSA can find out the strong sea-ice-loss area during warming hiatus, indicating a fundamental physical property change of the surface. It should be noted that the discussion of PDF change during the global warming slowdown period is incomplete in this study. We just chose several typical sites to illustrate the abundant information from the significant PDF change and verify the capability of the WSA algorithm. A more detailed discussion on scattered significant signals in the Indian Ocean and South America is another interesting topic. Meanwhile, we only applied the WSA in SAT in this study. Similarly, we can also apply the WSA to precipitation, wind speed, and pressure, which may also offer more insights into climate change. However, the current WSA method could only detect whether the distribution has changed or not, and the direction of change is not well depicted. For example, Fig. 1(a) can only give us a clue to discover PDF change signals, but cannot directly tell us whether most of the regions have more cold extremes or warm extremes. Therefore, how to measure the direction of PDF change is still a challenging topic and needs to be explored in the future. This section gives an introduction to optimal transport and Wasserstein distance(Wdistance). In practice, W-distance can quantify the distance and similarities between two probability distributions. Compared to other similarity metrics like Kullback-Leibler divergence, Wasserstein distance is rigorously defined and satisfies the metric axioms (Figalli & Glaudo, 2021). Wasserstein distance provides a powerful metric to quantify the similarities and discrepancies between probability distributions, and has been applied to climate science, such as model evaluation (Vissio et al., 2020), oceanographic data analysis (Hyun et al., 2022), and data assimilation (Tamang, 2022). In 1781, Gaspard Monge proposed the concept of optimal transport from one practical situation: if one uses soil to build fortifications, what is the cheapest way to transport the soil? Let  and  be two probability measures (distributions) on R d , this scenario leads to the following Monge formulation of optimal transport inf where c(x, T (x)) is the cost of transporting unit mass from x to T (x), and  = T #  means that for any (measurable Such T is said to be a transport map, and since T is a map, the mass at x could only be transported to one destination, which means Monge formulation does not allow splitting mass. In the 1940s, Leonid Kantorovich revisited Monge's problem and gave relaxation to Monge's formulation by splitting the mass. Kantorovich relaxation leads to the following formulation of optimal transport inf where   (, ) means  is a joint distribution with marginals  and , i.e.
When the cost function c(x, y) = |x  y| p , the metric side of Kantorovich formulation makes it valid to quantify the similarities between probability measures  and  via p-Wasserstein distance, defined as below where p  [1, ). In this work, we use the following earth mover's distance, i.e. p = 1, Earth mover's distance can be seen as the minimum amount of "work" required to transform mass from the probability measure (distribution)  into another probability measure (distribution) , where the "work" is measured as the amount of distribution weight that must be moved, multiplied by the distance it has to be moved. The code to calculate the earth mover's distance between two 1D distributions has been integrated into the wasserstein distance function of Scipy package. The input of wasserstein distance is just the two sequences of observed values in the empirical distributions, then the function returns the computed distance between these two distributions. The details can be found at (Scipy, 2023). Another option to compute Wasserstein distance is to use the ot.emd2 function in Python optimal transport package P OT (Flamary et al., 2021). The ot.emd2 function could compute the earth mover distance in 1D as long as the cost matrix is given by the corresponding l 1 distance. While P OT is one of the most efficient exact optimal transport solvers, it has not been designed to handle large-scale optimal transport problems. Therefore, if one needs to solve optimal transport with large number of samples, we do not recommend to use P OT . Interested readers may refer to P OT homepage for more information. Algorithm 1 provides details for the first step W-distance test in Wasserstein Stability Analysis. Given two samples X and Y , both series are first normalized individually by the following scalar x    where  is the mean value and  is the standard deviation. Then Wasserstein distance W 1 (X, Y ) between X and Y are obtained after normalization. It is worth noting that this normalization is necessary since the geographic difference may lead to a large variation in W-distance. For example, the high surface temperature variation at high latitudes usually leads to a larger W-distance than its counterpart in tropical zones. Nevertheless, it does not necessarily mean the PDF change at high latitudes is stronger than that at low latitudes. Furthermore, the length of the two samples X and Y in WAS could be different. Our significance test is based on the Monte Carro test (Xie et al., 2017). The null hypothesis H 0 of WSA is that the PDF discrepancy between two samples X and Y can be explained by the white noise. Clearly, the alternate hypothesis (H 1 ) is that the discrepancy of PDFs between X and Y results from factors other than white noise. After getting Wasserstein distance W 1 (X, Y ), we generate two new samples, X N and Y N , by adding white noise N (0, 1) to X and Y , where N (0, 1) is the standard Gaussian distribution. Then one could obtain the W-distance W 1 (X N , Y N ) between the new samples X N and Y N . Such additional operation of white noise is performed 500 times. Therefore, one could get a sequence of W-distances {W 1 (X N , Y N )} 500 N =1 , and then get a confidence interval (C.I.) of W-distances. Under the null hypothesis H 0 that the distribution discrepancy between X and Y can be explained by the white noise N (0, 1), the W-distance between the original samples X and Y must be indistinguishable from the counterpart distance between perturbed samples X N and Y N , i.e., W 1 (X, Y ) must be indistinguishable from W 1 (X N , Y N ). Under the null hypothesis H 0 with the confidence level 1  , one would expect that W 1 (X, Y ) is within the confidence interval of {W 1 (X N , Y N )} 500 N =1 . If not, H 0 is rejected, and one could claim that the distribution discrepancy between X and Y is not from the white noise N (0, 1) with the confidence level of 1  , where  = 0.01. In this study, a higher confidence level of 1   = 99%, rather than the common 95%, is used, since W-distance is very sensitive to the PDF change. Furthermore, when adding white noise to the normalized samples of X and Y , we use the standard Gaussian distribution N (0, 1). Since a large size of samples is used, according to the central limit theorem, the amplitude of white noise is set as one standard deviation individually. After performing the W-distance test for all locations, we could detect the places where the W-distance is significant, which are regarded as the significant PDF change zones. After getting the significant PDF change zones, we could do PDF analysis: plot two PDFs and the corresponding change for each significant change zone, to understand the climate change mechanism which makes W-distance significant. Algorithm 1 Wasserstein Stability Analysis (WSA) 1: procedure (H 0 : the discrepancy in distributions between X and Y is from white noise, with confidence level  = 0.01) 2: X and Y are normalized individually by the x  scalar.
Get W 1 (X, Y ) 4: N  1 5: while N  500 do 6: X N  X + random(0, 1)  add white noise N (0, 1) to X and Y Y N  Y + random(0, 1) Get W 1 (X N , Y N )  X N and Y N are series with white noise 9: N  N + 1 10: end while reject H 0  The discrepancy is not from N (0, 1) 13: else 14: H 0 is not rejected 15: end if 16: end procedure The authors greatly appreciate the comments of Dr. Zhaolu Hou from the Ocean University of China and Dr. Shuailei Yao from the Institute of of Sciences. This work is The authors greatly appreciate the comments of Dr. Zhaolu Hou from the Ocean University of China and Dr. Shuailei Yao from the Institute of of Sciences. This work is Surface 2 meters air temperature data is from ERA5 reanalysis products of the European Centre for Medium-Range Weather Forecast (ECMWF) (Hersbach et al., 2020), which can be accessed at the following link Sea ice data is from the and Sea Surface Temperature dataset (HadISST) of Met Office in British (Rayner et al., 2003), which could be accessed via the following link The code of Wasserstein Stability Analysis can be found online via the link This material is the supporting information for Wasserstein distance in optimal transport and W-distance test in Wasserstein Stability Analysis. In the first section, We give an introduction to Wasserstein distance, and especially we list the open source code to compute Wasserstein distance. Then in the second section, we give a detailed description of the W-distance test of Wasserstein Stability Analysis, which is formulated in Algorithm 1 Surface 2 meters air temperature data is from ERA5 reanalysis products of the European Centre for Medium-Range Weather Forecast (ECMWF) (Hersbach et al., 2020), which can be accessed at the following link Sea ice data is from the and Sea Surface Temperature dataset (HadISST) of Met Office in British (Rayner et al., 2003), which could be accessed via the following link The code of Wasserstein Stability Analysis can be found online via the link This material is the supporting information for Wasserstein distance in optimal transport and W-distance test in Wasserstein Stability Analysis. In the first section, We give an introduction to Wasserstein distance, and especially we list the open source code to compute Wasserstein distance. Then in the second section, we give a detailed description of the W-distance test of Wasserstein Stability Analysis, which is formulated in Algorithm 1