A wavelet transform method to determine monsoon onset and retreat from precipitation time‐series

A new method to determine monsoon onset and retreat timings using wavelet transform methodology applied to precipitation time‐series at the pentad scale is described. The principal advantage of this method is its portability, since it can be easily adapted for any region and dataset. The application of the method is illustrated for the North American Monsoon and the Indian Monsoon using four different precipitation datasets and climate model output. The method is shown to be robust across all the datasets and both monsoon regions. The mean onset and retreat dates agree well with previous methods. Spatial distributions of the precipitation and circulation anomalies identified around the onset and retreat dates are also consistent with previous work and illustrate that this method may be used at the grid‐box scale, not just over large area‐averaged regions. The method is also used to characterize the strength and timing of the Midsummer drought (MSD) in southern Mexico and Central America. A two peak structure is found to be a robust structure in only in 33% of the years, with other years showing only one peak or no signs of a bimodal distribution. The two‐peak structure analysed at the grid‐box scale is shown to be a significant signal in several regions of Central America and southern Mexico. The methodology is also applied to climate model output from the Met Office Hadley Centre UKESM1 and HadGEM3 CMIP6 experiments. The modelled onset and retreat dates agree well with observations in the North American Monsoon but not in the Indian Monsoon. The start and end of the modelled MSD in southern Mexico and Central America is delayed by one pentad and has a stronger bimodal signal than observed.


| INTRODUCTION
Monsoon regions are characterized by a strong contrast between the rainy and the dry seasons, with at least 70% of the total annual rainfall observed during summer in a monsoon (Zhou et al., 2016;Wang et al., 2017). The start and end of the wet season, defined by onset and retreat (or demise) dates, are physically associated with sudden changes in the moisture transport driven by the changes to the circulation. Around the onset and retreat dates, changes in the predominant wind direction or strength also cause abrupt changes in meteorological conditions such as rainfall, humidity and cloud cover (Zhou et al., 2016;Gadgil, 2018). A current view is that monsoons can be defined based solely on the contribution of summer rainfall to the annual total highlighting the importance of rainfall to several societal sectors (Zhang and Wang, 2008;Wang et al., 2017).
The timing and strength of the rainy season are key aspects of the seasonal cycle of monsoon regions such that the length of the rainy season greatly influences sectors such as agriculture (Sultan et al., 2005;Gadgil and Rupa Kumar, 2006;Jain et al., 2015;Harvey et al., 2018) and water management (Turner and Annamalai, 2012;Bussmann et al., 2016). The objective determination of onset and retreat dates is key for climate and weather research aiming to understand the short and long-term variability, trends and predictability of the regional monsoons (e.g., Kitoh and Uchiyama, 2006;Cook and Buckley, 2009;Htway and Matsumoto, 2011;Lucas-Picher et al., 2011;Nieto-Ferreira and Rickenbach, 2011).
For this reason, a plethora of methods have been used to diagnose the onset and retreat dates from a range of variables and datasets. Bombardi et al. (2020) provides a recent review of methods for objectively defining monsoon onset and retreat dates and highlighted the technical differences and purposes of each. Methods can be divided into those that aim to evaluate monsoon onset and retreat on a regional scale (e.g., Webster and Yang, 1992;Fasullo and Webster, 2003;Garcia and Kayano, 2013) or at a local or grid-box scale (e.g., Liebmann and Marengo, 2001;Cook and Buckley, 2009).
Threshold methods are the most commonly used local-scale methods that typically diagnose onset and retreat from a precipitation time-series (Bombardi et al., 2020). These methods evaluate the accumulated  or daily/pentad-mean rainfall rates (Geil et al., 2013) and determine the onset and retreat dates when the time-series exceeds or falls below a pre-defined value (threshold) for a given amount of time (persistence). The persistence parameter is used to reduce the method's sensitivity to noise in the precipitation time-series. The threshold parameter can be a statistical measure of the seasonal cycle such as the total annual mean rainfall (Arias et al., 2012) but the threshold and persistence parameters can also be subjectively determined. The parameters for threshold methods vary distinctly from study to study for several reasons. First, the parameters vary from one monsoon region to another because the wet and dry seasons in each region have different timings, strengths and dynamical features (Wang et al., 2017). For instance, the North American Monsoon occurs over a much shorter time than the South American Monsoon and is smaller geographically (Vera et al., 2006). Within a given monsoon region, for example, in the South American Monsoon, several methods are used for different purposes depending on the temporal and spatial scales of interest (see, for example, Liebmann and Marengo, 2001;Marengo et al., 2001;Nieto-Ferreira and Rickenbach, 2011;Carvalho et al., 2011;Garcia and Kayano, 2013). The differences between the characteristics of the datasets, for example, in horizontal resolution, relate to differences in the seasonal cycle of rainfall which means that the implementation of threshold methods in different datasets also requires normalization or statistical treatment of the threshold and persistence parameters.
In other words, each threshold method is tailored to a monsoon region using a specific dataset and a specific variable for a given purpose. This characteristic of the threshold methods poses various shortcomings. Firstly, practical shortcomings of the threshold methods, particularly rigid thresholds, include false hits (Moron and Robertson, 2014) or some years not meeting the threshold and persistence criteria (Arias et al., 2012) requiring further relaxation of the parameters. Second, the conclusions of most analyses of threshold methods arise from the output of one single-method and one single dataset. These studies would benefit from the input of a second independent method to confirm the results. However, due to the differences between regional monsoons most techniques are not suitable to be used in other regions.
The Coupled Model Intercomparison Project (CMIP) assessments of monsoon onset and retreat typically use precipitation threshold methods due to the lack of high temporal or vertical resolution output from all models to estimate vertically integrated quantities need for some methods (e.g., Geil et al., 2013;Zou and Zhou, 2015;Ha et al., 2020). Threshold methods have multiple shortcomings for CMIP assessments as the persistence and threshold parameters are tuned for observations with a specific seasonal cycle but models have a range of biases in the seasonality, magnitude and spatial distribution of rainfall (Pascale et al., 2019;García-Franco et al., 2020). The use of pre-defined threshold values may also not be suitable to compare different model experiments with changes in forcing where the climatological mean rainfall or the seasonal cycle may change within the model run. These shortcomings are relevant because a proper diagnosis of the seasonal cycle in CMIP assessments is key to understand and diagnose current and future changes to monsoon seasonality as a result of greenhouse forcing (Zhou et al., 2016;Wang et al., 2017). The observational or model-based analysis of the monsoons as a global phenomena also require a systematic approach to determine onset and retreat dates that is not tailored to any specific region (Zeng and Lu, 2004). For example, projects such as the Global Monsoons Model Inter-comparison Project (GMMIP) which is part of CMIP6 may benefit from broader and more systematic methods to be used across models and monsoon regions (Zhou et al., 2016) and particularly methods that can be implemented with standard output from modelling centres. These systematic or large-scale methods may also be relevant in the evaluation of trends in the timings of the rainy season on hemispheric or global scales in observations and models (Zeng and Lu, 2004;Zhang and Wang, 2008).
The objective diagnosis of shorter time-scale rainfall variability, such as bimodal regimes and active and break phases of a monsoon, also requires methods that can separate relatively drier and wetter periods within the rainy season. For instance, the characteristics of the seasonal cycle of precipitation in northwestern Central America and southern Mexico fit the definition of a monsoon climate (Wang et al., 2017), characterized by the majority of precipitation observed during local summer, and for this reason this region is sometimes considered to be part of the American monsoon (Vera et al., 2006;Pascale et al., 2019). However, this region shows a unique climatological bimodal signal (Karnauskas et al., 2013). After monsoon onset, rainfall decreases considerably around the midsummer; this decrease is followed by a secondary increase in precipitation in the late summer (Mosiño and García, 1966;Magaña et al., 1999).
Several studies have documented the characteristics of this bimodal signal, most commonly referred to as the Midsummer drought (MSD; e.g., Mosiño and García, 1966;Karnauskas et al., 2013;Perdig on-Morales et al., 2018) and have linked the cause of the bimodal distribution of rainfall to regional changes in the circulation and the moisture transport (Magaña et al., 1999;Herrera et al., 2015). However, the objective determination of the strength, spatial distribution and robustness of bimodal signals is not straightforward. For example, the global method used in Bombardi et al. (2020) fails to diagnose the region of southern Mexico, Central America and the Caribbean as a bi-modal regime.
The majority of existing methods to diagnose bimodal signals in this region use geometric or statistical measures of the monthly mean rainfall that measure the difference between the months of maximum rainfall and the drier months. However, this approach fails to capture the shorter-scale changes that have been shown to occur in both observations and model data, as the MSD does not start or end exactly on given calendar months (Magaña et al., 1999;García-Franco et al., 2020). Zhao and Zhang (2020) review and compare several methods to detect and measure the MSD, finding that using monthly-mean data and prior assumptions of the dates of the first and second peaks can lead to errors. Only a handful of methods exist that can determine the characteristics of the MSD in Central America on sub-monthly timescales. Anderson et al. (2019) analysed the pentadmean time-series from the Climate Hazards Infrared Precipitation with Stations (CHIRPS) dataset. After a double temporal smoothing of the time-series, Anderson et al. (2019) determined the MSD timings through a threshold method with parameters tailored to the CHIRPS dataset. In turn,  used daily mean time-series and determined the two-peak structure through linear-regression analysis.
In short, research on monsoon regions requires objective methods to determine the timings of onset, retreat and short-scale changes within the rainy season such as the MSD and active and break phases. Multiple methods exist, each with various parameters fit for different purposes, but these methods present shortcomings for studies that compare results from multiple datasets or investigate model experiments where climatological rainfall and the seasonal cycle change. Similarly, studies that investigate the impact of decadal modes of variability (e.g., Arias et al., 2012), greenhouse warming (e.g., Geil et al., 2013) or general trends (e.g., Sahana et al., 2015) rely solely on the output of one single method whereas the use of two or more methods may help to test the sensitivity of their results to the chosen method and parameters. Both the objective determination of monsoon onset and retreat and the timings of bimodal regimes require a method that can analyse temporal changes to precipitation on several scales and that can be used on any gridded dataset.
The goal of this study is then to present an objective approach that is more portable across datasets, regions and robust for various purposes. This article introduces a wavelet transform method to determine monsoon onset and retreat dates using pentad-mean precipitation timeseries. Wavelet algorithms have been extensively used in atmospheric research for multiple purposes, such as the detection of the boundary layer height (e.g., Brooks, 2003), as well as to analyse time-frequency features of a signal (e.g., Whitcher et al., 2000;Dimdore-Miles et al., 2021). In fact, Allen and Mapes (2017) used a wavelet analysis to determine monsoon onset and retreat using daily OLR data. The method is constructed such that the determination of monsoon onset and withdrawal dates is less sensitive to the characteristics of the time-series, that is, the characteristics of the seasonality of each monsoon or of a given observational dataset. Furthermore, the method is expanded to characterize bimodal regimes which is illustrated for the MSD of southern Mexico, Central America and the Caribbean.
The remainder of this article is organized as follows: Section 2 describes the methods and datasets. Section 3 shows the results of applying the method to the Indian and North American Monsoons and the MSD. Section 4 summarizes the method and discusses the results.
2 | DATA AND METHODS 2.1 | Data 2.1.1 | Observations and reanalysis data Table 1 summarizes the datasets and variables used in this study. We use three gridded precipitation datasets, the Tropical Rainfall Measurement Mission (TRMM) v7 3B42, Climate Hazards Infrared Precipitation with Stations (CHIRPS) and the Climate Prediction Center Merged Analysis of Precipitation (CMAP). These three precipitation datasets are merged products, TRMM and CMAP mainly use microwave satellite measurements complimented by several other sensors and calibrated with rain-gauge data whereas CHIRPS uses several products from TRMM, as well as high-resolution station data. These datasets also differ in their end-product horizontal resolutions.
The precipitation output from the latest reanalysis from the European Centre for Medium-Range Forecasting: ERA-5, which has been shown to exhibit a relatively good representation of the temporal characteristics of rainfall in monsoon regions (e.g., García-Franco et al., 2020). Other variables from ERA5 used to diagnose changes to the circulation associated with monsoon onset were geopotential height at 500 hPa and wind speed (u ! ).

| Model data
Daily precipitation data from the CMIP6 archive are used, retrieved from: https://esgf-index1.ceda.ac.uk/ projects/cmip6-ceda/, to illustrate the method using standard climate model output. In particular, we use results from the Met Office Hadley Centre models (MOHC): HadGEM3 GC3.1 and UKESM1 and we have chosen the pre-industrial control and historical experiments Sellar et al., 2019;Andrews et al., 2020) which are amongst the most commonly used experiments in CMIP6. The daily precipitation data were converted to pentad-scales. Table S1 provides a summary of the experiments used and the acronyms of each experiment for each model Ridley et al., 2018Ridley et al., , 2019Sellar et al., 2019). HadGEM3 GC3.1 was run at two horizontal resolutions, N96 and N216 for the pre-industrial control experiment (Table S1). UKESM1 differs from HadGEM3 by representing Earth System processes, such as interactive ocean biogeochemistry and atmospheric chemical interactions (Sellar et al., 2019). The pre-industrial Control experiments use constant forcing estimated for 1850 whereas the historical experiments aim to represent timevarying greenhouse emissions, volcanic eruption and solar signals in the historical period (1850-2014) Eyring et al. (2016). These simulations were recently evaluated in the North American Monsoon and the MSD region, showing a good representation of the seasonal cycle in these regions (García-Franco et al., 2020).

| Wavelet transform method
Wavelets are band-limited wave-like functions with specific mathematical properties that include a finite energy and zero-mean (Whitcher et al., 2000;Addison, 2017).  (2017) Note: For each dataset, the acronym used hereafter, the period of coverage, the field used and the horizontal resolution are shown.
The wavelet function is defined using two parameters, a dilation (a width or temporal scale) and a translation (centroid in time/space). Wavelet transforms are the result of the inner product (convolution) of a wavelet function with a time-series or a signal (Addison, 2017). The wavelet transform can be thought of as a local comparison between the wavelet function and the observed signal for different frequencies or temporal scales. The information provided by a wavelet transform largely depends on the characteristics of the wavelet function used, thus, different wavelet functions are used for different purposes (Addison, 2017). For the purpose of finding the onset and retreat dates, the wavelet based on the Haar function is useful as this wavelet is suitable to find sudden changes in a signal (Addison, 2017). The Haar wavelet is defined as the noncontinous piece-wise function: where, a is the dilation coefficient, b is the centre of the wavelet or the translation coefficient and t is the time coordinate.
The wavelet covariant transform is then the inner product of the Haar wavelet with a timeseries (Brooks, 2003), that is: where, pr(t) is a time-series of precipitation, either on daily or pentad scales and W f (a, b) is the matrix of the covariant transform and t i and t f are the start and end time-points. No statistical treatment, normalization or anomaly, a priori, is calculated on the precipitation timeseries pr(t) so the units of W f are the same as the precipitation time-series (e.g., mm day −1 ). Monsoon timings can be observed as sharp changes to precipitation, that is, rainfall sharply increases at onset and sharply decreases at retreat. However, measuring these changes can be difficult since precipitation timeseries are typically noisy. The Haar wavelet is useful in these cases for signal detection since the wavelet transform is interpreted as gradients across different temporal scales and can smooth out the high-frequency variability using sufficiently large dilation scales. In other words, the wavelet covariant transform (W f (a, b)) measures gradients on a dilation scale a for each time-step (b). Figure 1 shows the Haar wavelet and 1 year of observed precipitation in the North American Monsoon from the CMAP dataset. Figure 1(a) illustrates how the wavelet function compares the observed signal in the interval b<t≤b + a 2 with the values of the signal in the interval b − a 2 ≤t<b where b in this case is a pentad time step. The wavelet transform coefficient for dilation a = 20 pentads at the translation of b = 35, that is, pentad 35, is a measure of the precipitation difference between the sum of the observed rainfall 10 pentads after pentad 35 and the sum of the observed rainfall 10 pentads before pentad 35 as illustrated in Figure 1(b). Figure 2 shows an example of the wavelet transform (WT) application using the observed climatological precipitation in the North American Monsoon in four different precipitation datasets. The mean climatological rainfall rates (upper panel) differ in their peak summer rainfall rates but qualitatively show similarities in the start and end dates of the rainy season. The wavelet transform coefficients (W f (a, b) in the middle panel) for a small dilation a = 4 are relatively noisy but show a clear maximum and minimum that correspond well with the maximum and minimum of longer dilations (a = 14, 20). The sum of these four coefficients at each translation or pentad b, highlight a maximum found around June 22 and a minimum found around September 21st, which agree well with previous results of mean onset and retreat dates in the North American Monsoon (e.g., Arias et al., 2012;Geil et al., 2013).
Local maxima in the wavelet transform highlight positive steps in the precipitation time-series with a coherent scale of a pentad steps. This interpretation is then extended to diagnose monsoon onset. The pentad (b*) corresponding to the maximum of the sum of the transform over a set of scales is defined as monsoon onset (MO), that is: where, a 0 and a f are the limits of the pentad scales, that is, the dilation coefficients, b * is the pentad of maximum W f (a, b) P W f (a, b) and the monsoon onset pentad. Similarly, the monsoon retreat pentad is found at the minimum of the sum of the wavelet transforms, that is, In other words, we seek to find monsoon onset and retreat using the maximum and minimum the wavelet power spectrum over a range of temporal scales. Several sensitivity tests were performed with different dilation coefficients (a) in the different observational datasets, models and regions and a set or vector of dilation scales was found to be optimal to be used for all purposes. The set of dilation coefficients a ! = 28,30, ÁÁÁ,54 ð Þwas found to be robust, that is, was able to capture the onset and retreat dates in all the datasets.

| Identification of monsoon onset and retreat
Monsoon onset is defined as the maximum sum of wavelet coefficients, capturing positive gradients within the scales of 14-27 pentads (half of the elements of vector a defined above). Monsoon Retreat has a similar definition, capturing the greatest negative gradient of precipitation over the same pentad scales.
For example, Figure 3 illustrates the method in the North American Monsoon in the TRMM dataset for 2009. Figure 3a shows the WT coefficient matrix, showing the changes in precipitation for dilations ranging from 10 to 40. A clear signal of positive coefficients is observed between pentads 28-34 and a similar negative signal observed in pentads 56-60. Figure 3b shows the timeseries of the observed precipitation, which suggests that monsoon onset occurs sometime between pentads 28 and 34. Observed rainfall rapidly decreases after pentad 59 suggesting that monsoon retreat can be diagnosed around this pentad. Figure 3c shows the WT coefficients as a function of pentad for several dilations (a 0 ). The coefficients for all scales seem to follow a very similar behaviour, increasing during spring to reach a maximum around pentad 32 and thereafter decreasing to a minimum around pentad 59. When the sum of the wavelet transform coefficients across the dilations is computed (Figure 3d) this behaviour becomes much clearer. The maximum and minimum are found at pentads 32 and 59, respectively and these pentads define the onset and retreat times. For comparison, the results from the method of Geil et al. (2013) are shown in Figure 3b, indicating that this method may have found an earlier onset. Figure S1 shows another example in the same region but using model data from HadGEM3 GC3.1N216.

| Extension for application to the MSD signal
The wavelet method can be extended to characterize the shorter scale variations of precipitation of the MSD in Central America and the Caribbean. First, the monsoon onset and retreat dates are determined in the time-series from the area-averaged precipitation in the MSD region via the approach described in the previous section. Once the onset and retreat dates are established, an additional wavelet analysis determines the dates in which the MSD starts and ends. The onset and end of the drier period of the MSD can be found by computing the wavelet transform again but using smaller dilations and over a limited temporal range. In particular, the WT is only calculated in the 20 pentads before and after the dates defining monsoon retreat and onset, respectively. The MSD Onset (MSDO) and MSD End (MSDE) are defined as the minimum and maximum, respectively, of the sum of the wavelet transforms (Equations 3 and 4) using dilation coefficients a ! Ã = 10,12, ÁÁÁ,24 ð Þ . Figure 4 illustrates the use of the WT method to determine the dates of MSDO and MSDE for the precipitation of 2017 in ERA-5. Figure 4a depicts the wavelet covariant transform matrix, showing the W f coefficients for each dilation a at each pentad b. The onset of rainfall is diagnosed around the time-steps of highest positive W f coefficients-around pentads 24-32 for almost all dilations. These positive coefficients are followed by a period of negative values from pentad 32 to pentad 40, which represent the decrease in precipitation, or relative drought, in the midsummer. The MSD is followed by another period of positive coefficients from pentad 44 to pentad 52, illustrating the so-called second peak of precipitation and, finally, a period of negative coefficients associated with monsoon retreat.
The coefficients of the wavelet transform [W f (a 0 , b)] for selected dilations a 0 (Figure 4c) show that the smaller dilations are more sensitive to smaller scale variations in the time series and longer dilations better highlight the long-term change of the time series. For example a 0 = 18 shows signs of a MSD by showing two local maxima and two local minima, whereas a 0 = 54 only shows a local maximum and a local minimum associated with onset and retreat.
The maximum and minimum of the sum over all dilations (Figure 4) depict the rainfall onset and retreat dates, respectively. The second wavelet transform W f (a * , b) is computed over smaller dilation coefficients (a * ) near the onset and retreat dates as described above to highlight the MSDO and the MSDE. Figure 4e shows the sum of the wavelet transform coefficients W f (a * , b) and the pentad of the MSDO, 34, and MSDE, 49, corresponding to the minimum and maximum of the sum of these wavelet transform coefficients, respectively.
The strength of the MSD can be measured through the maximum and minimum sum of the coefficients P W f (a * , b) used to define the start and end dates of the MSD. For example, in Figure 4(e) the minimum of the P W f (a * , b) was −20 mm day −1 found at pentad 35 and an opposite local maximum of +20 mm day −1 at pentad 49. These two values, hereafter coef1 and coef2, provide a quantitative measure of the strength of the MSD for this year in this dataset and will be used to measure the spatial variability of the magnitude of the MSD in the different datasets.

| Comparative methodologies
For validation purposes the wavelet transform method is compared with existing methods that determine onset and retreat in the North American and Indian monsoons. The threshold methods of Geil et al. (2013;hereafter G13) and Arias et al. (2012;hereafter A12) are compared with the results of the wavelet transform in the North American Monsoon. G13 used a threshold of 1.3 mm day −1 for at least 3 days for onset and 7 days for retreat for daily TRMM observations. In this study we adapt this method for TRMM to the same threshold value, but the onset pentad is the first pentad above the threshold whereas for retreat, we require rainfall to be below the threshold for two consecutive pentads. The method of A12 defines onset with two conditions. The first condition to find the onset pentad is that six out of the eight subsequent pentads must have rain-rates above the annual-mean climatological rainfall. The second condition is that at least six out of the eight previous pentads must be below the annual-mean climatological rainfall. The opposite definition is used to determine the pentad of monsoon retreat.
In the Indian Monsoon region, a commonly used method is known as the hydrologic onset and withdrawal index (HOWI) based on moisture transport over the Arabian Sea (Fasullo and Webster, 2003;Sahana et al., 2015;Chevuturi et al., 2019). To compute the HOWI index, first, the vertically integrated moisture transport (VIMT) is computed from daily ERA-5 data in the Arabian sea, as described by Fasullo and Webster (2003), that is: where, χ is the VIMT, g is the gravitational acceleration, p are the pressure levels, q is the specific humidity and V is the wind vector. The VIMT is then normalized using the transformation: where, χ is the unnormalized time-series, χ is the mean seasonal cycle of the unnormalized index and HOWI is the normalized index. The onset date is defined as the first day of each year where the HOWI index is greater than zero and the retreat date is the first day after the onset date that the HOWI index is negative (Fasullo and Webster, 2003;Sahana et al., 2015). The necessary daily data of moisture and wind speed on sufficient vertical levels to compute the HOWI index in the MOHC submissions to CMIP6 was not available, so the HOWI index can only be computed using ERA5 and will be compared with the WT method used on the observational gridded datasets.

| RESULTS
The onset and retreat dates were determined for each year in each observed and model dataset for the Indian, North American and MSD regions using the methods described in the previous section. The calculations were performed for area-averaged precipitation time-series representative of the core regions defining these monsoons. Calculations were also made at grid-box scales to illustrate the spatial distribution of the onset and retreat dates. Table 2 shows the onset and retreat dates of the North American Monsoon estimated using the G13, A12 and WT methods. The table reports the results for three observational datasets, ERA-5 reanalysis and five climate model experiments. The observations generally agree that the onset date is found at pentad 33 (around June 15), according to the WT and the method by G13. However, the method of G13 generally places the retreat date one pentad after the WT method, that is, around October 7th for G13's method and October 2nd for the WT's method. The method by A12 disagrees with G13 and the WT methods on both onset and retreat mean pentads, in both cases finding an earlier onset (pentad 30) and retreat (pentad 54). The climate models reasonably represent the mean onset and retreat dates, as only the onset dates from two experiments of UKESM1 are statistically different from two of the observational datasets. The similiarities in onset and retreat dates confirm that the seasonal cycle is very well represented by these models in the North American Monsoon. The method by A12 in the simulations also produces an earlier onset and retreat dates when compared with the other two methods by about 1.5 pentads, but this is within the uncertainty range given by the interannual variability. Figure 5 compares the estimated onset and retreat dates using these three methods in the North American Monsoon in 2010 using the three observational datasets and ERA-5. The method by A12 shows an earlier onset and retreat in all the datasets whereas the WT and G13 agree in almost every dataset. The WT method is the only method that estimated the same onset date for all the four datasets  and the same retreat date in three of the four datasets, with the fourth dataset delayed by only one pentad. The meteorological changes associated with onset and retreat in the North American Monsoon are shown in Figure 6. The composite differences of the precipitation, wind and geopotential changes 10 days prior and after onset and retreat are compared for the three methods. The impact of monsoon onset in precipitation is slightly stronger using the method by A12 than the WT or G13. The WT method shows a very similar pattern and magnitudes of the anomalies of onset and retreat when compared with the other two methods. The method by G13 produces the weakest anomalies, particularly of precipitation, but very similar patterns overall. Figure 7 shows the spatial distribution of the mean onset and retreat dates in the North American Monsoon region for various datasets. There is high agreement between TRMM, CHIRPS and ERA5 on the spatial pattern of mean onset and retreat dates. Onset in western Mexico is around pentad 31 (around June 1st) whereas in Chihuahua and Sonora, the rainy season begins shortly after pentad 35 (June 22nd). The pattern in the mediumresolution simulation GC3 N216 piControl is consistent with observations, particularly during onset. However, the spatial pattern of the mean retreat dates in the northern regions of the monsoon show an earlier than observed retreat, possibly associated with the dry bias in this region in these models (García-Franco et al., 2020). Table 3 compares the mean onset and retreat dates in the Indian Monsoon computed from the HOWI index using ERA5 data and the WT used for gridded precipitation datasets. The onset and retreat dates from the HOWI index were converted from the daily to the pentad-scale to compare with the WT. The mean onset date for the HOWI index is May 27th between pentads 29 and 30, and retreat is between pentads 49 and 50, around September 3rd. The mean onset date found using the WT method for the four observational datasets was close to pentad 32, about two pentads later than the HOWI index. The mean retreat date for the WT method was about one pentad earlier than the HOWI results. Overall, the models exhibited later than observed onset and retreat dates (4 pentads). The differences between the hydrological determination of onset and retreat dates, through HOWI, and the WT method on gridded datasets is significant, according to a Welch's t-test done between HOWI and all the gridded datasets. These differences may be due to the different regions where each method is defined, that is, HOWI is defined over the whole of the Arabian Sea where an earlier onset and later retreat would be expected when compared with rainfall over mainland India, where the WT method was applied. Figure 8 shows differences between the WT (based on ERA5 precipitation) and HOWI methods, comparing precipitation and moisture fluxes at 850 hPa 10 days prior to and following monsoon onset and retreat. The HOWI index better captures the moisture transport in the Arabian Sea whereas the WT method best captures precipitation differences prior to and following monsoon onset and retreat over mainland India. The WT method is also able to capture onset and retreat dates and the associated anomalies within the climate model output. Figure S2 shows the precipitation and moisture transport anomalies around the onset and retreat of the Indian monsoon in three different climate model experiments. While the models show significant biases in the timings of the monsoon, according to a Welch's t-test (Table 3) the models show similar patterns of rainfall and moisture transport anomalies to the observations around both onset and retreat.

| The Indian monsoon
In the case of the Indian monsoon, Figure 9 shows that the mean onset and retreat dates vary greatly spatially on the southern tip of the subcontinent. While most of northern India has a mean onset date around pentad 33, the western coast shows an earlier onset by about one or two pentads. There is high agreement in the onset date in the three observed/reanalysed datasets over mainland India and between TRMM and ERA5 over the western coast of India. Earliest onset is found on the western coast around pentad 25-27 and extending to central India by pentad 31. The GC3 N216 simulation, however, shows a later onset than observed by about two pentads in most regions. In contrast, the spatial pattern for the mean retreat date shows higher spatial

| The midsummer drought
Results from the application of the wavelet transform to the MSD, including the mean onset and retreat dates as well as the start and end dates of the MSD period, are reported in Table 4. The mean onset date in the observations is around pentad 27 (May 14th), whereas the retreat date is around pentad 61 (October 31st). The end of the so-called first-peak period, or start of the relatively drier period (MSD), is consistently found in all the observed datasets to be around pentad 36 (around June 29th). The end of the drier period or start of the second peak is also consistently determined to be between pentads 48 and 49 in the four observational datasets. In other words, the MSD has a mean duration of 12 pentads, or around 2 months, from late June to late August. In the MOHC simulations, the MSD starts slightly later than observed by about two pentads, and ends about one pentad later than observed. Figure 10 shows the rainfall anomalies associated with the different periods (stages) of the rainy season in southern Mexico and Central America. These include monsoon onset and retreat, and the start and end of the MSD, the MSDO and MSDE, respectively. For each stage, we compared the anomalies computed by separating the stages using the WT method or the dates of the climatological monsoon onset, retreat, MSDO and MSDE as found in Table 4. In this way, the ability of the WT method to characterize rainfall variations is tested against a first best guess-the climatological mean dates. Overall, using the dates for MSDO and MSDE from the climatological dates results in weaker anomalies than compositing via the specific dates for each year obtained with the WT method. Even though the area-averaged signal used diagnose the different MSD stages focuses on a small region of southern Mexico and northern Central America, the anomalies associated with the onset and end of the MSD (Figure 10f,h) extend across the East Pacific warm pool, most of the western coast of Mexico and into to the Caribbean Sea and Cuba. The analysis of individual years of observed precipitation in the selected area-averaged time-series showed that not all years showed a bimodal signal in the area-averaged precipitation (see Figure S3). In fact, a given year could be classified as having (a) a canonical two-peak structure separated by an MSD, (b) only having a first peak and an MSD but no second peak, (c) only having a second peak but no clear MSD or (d) a plateau-like monsoon season with no MSD-type variations (see Figure S3).
Due to this year-to-year variability in the characteristics of the seasonal cycle, an objective measure was defined to determine whether a signal presented a robust MSD-bimodal signature. For this purpose, the WT algorithm was applied to randomly generated precipitation time-series.
The random time-series are constructed by randomly sampling observations in the wet and dry seasons. The pentad-mean onset and retreat dates from Table 4 were used to composite the observations into dry and wet distributions. For a random time-series, the values of each pentad are randomly selected from the dry or wet distributions, depending on the pentad. In this way, the value of a given pentad of the random time-series may have been observed at a different pentad; the only constraint is that the random values come from pentads that were observed in the same season. This approach has two advantages. First, that the approach imposes the monsoon-like feature of a sharp wet-dry contrast but secondly, the random selection in the wet season removes the possible signal of the MSD in the climatological rainfall.
The random time series are then constructed by randomly drawing values at each pentad from the wet or dry season distributions of each dataset. Then, the WT method was used on 10,000 of these random-time series. This approach rendered a distribution of coefficients (coef1 and coef2) essentially representing the variability of the WT method applied to noise. Figure 11 shows the pentad-mean time-series from 2 years in the TRMM dataset, and two randomly  generated time-series. The coefficients coef1 and coef2, illustrated in Figure 11b, measure the difference in precipitation between the first peak and the MSD period and the MSD and the second peak, respectively. The first quantile of coef1 and the third quantile of coef2 provide a measure of robustness for the observed coef1 and coef2.
In other words, for a year to be classified as having a robust MSD signal, the resulting coef1 and coef2 of the WT procedure must be lower and higher, respectively, than those obtained for a random time-series. The analysis of coef1 then determines the existence of a first-peak MSD type variability and coef2 determines the robustness of a possible second-peak for that year. By this procedure, a given year could fit into four categories: • Canonical MSD: coef1 lower than the first quantile (25%) of random coef1 and coef2 higher than the third quantile (75%) of random coef2. • First peak + MSD: coef1 lower than the first quantile of random coef1 but coef2 lower than the third quantile of random coef2. In other words, the second peak is not distinguishable from noise. • Second peak only: coef1 higher than the first quantile of random coef1 but coef2 higher than the third quantile of random coef2. In other words, the second peak is distinguishable from noise, but there is no first-peak + MSD structure. • No MSD: coef1 higher than the first quantile of random coef1 and coef2 lower than the third quantile of random coef2. In other words, the precipitation time-series shows no robust signal of an MSD regime, with a first or second peak. Figure S3 shows how separating years into these categories affects the pentad-mean seasonal cycle of precipitation in southern Mexico and Central America in four observational datasets. This figure also validates the above procedure as the WT method is able to robustly separate years into the different categories.
For each dataset we determine those grid-points showing a robust MSD. We use the climatological rainy and wet seasons as above, to construct the random timeseries for each grid-point and estimate the random values of coef1 and coef2, repeating the procedure 10,000 times. A given grid-point is diagnosed to have a MSD when the value of coef2 − coef1 is higher than the third quantile of the PDF of the random time series. The value of coef2 − coef1 is a measure of the magnitude of the MSD since coef2 measures the relative strength of the second-peak compared with the MSD and therefore positive in an MSD grid-point and coef1 compares how dry the MSD is relative to the first-peak and thus negative if an MSD regime is observed at that grid point. Figure 12 shows the regions where the climatological rainfall shows a MSD signal that is distinguishable from noise, that is, the values of coef2 − coef1 exceed the third quantile of the distribution composited with random timeseries. Cuba, western Central America and most of southern and central-eastern Mexico exhibit a robust MSD signal., This map also shows the MSD amplitude; the strongest MSD signal is found on the western coast of northern Central America and northeastern Mexico. The high correspondence between the three observational datasets shows that the method is robust across datasets. These results In particular, the method is able to replicate the previously reported MSD signal in the Pacific Mexican coast and the stronger MSD signal in northeastern Mexico. Figure 13 shows the spatial distribution of the mean onset and retreat pentads and the start and end of the MSD, in the grid-points where the signal is significant as in Figure 12. The earliest rainfall onset is found on the western coast of southern Mexico, Guatemala and El Salvador, as well as in Cuba, at pentad 25, whereas onset in the Yucatan peninsula is found at pentad 28 and even later, around pentad 31, in the eastern states of Mexico. In contrast, the retreat date seems spatially more homogeneous as northern Central American has a mean retreat date around pentad 59 and central Mexico around pentad 54. The MSD coherently starts over the western coast of Guatemala and Chiapas around pentad 33. In contrast, the MSD on the eastern Mexican states of Veracruz and Campeche begins after pentad 40. The earliest MSD end (Figure 13h-j) is found in central and northeastern Mexico, around pentad 42 whereas the MSD in Guatemala ends around pentad 48.

| SUMMARY AND DISCUSSION
A novel method is presented to diagnose monsoon onset and retreat dates using pentad-mean precipitation based on the computation of a wavelet transform over multiple temporal scales. The wavelet function used is the Haar wavelet, a wavelet typically used to find abrupt changes in signals. Onset is defined as the maximum of the sum of the coefficients of the wavelet transform computed over a range of temporal scales or dilations. These dilations were found to provide the best results in a range from 28 to 54 pentads. Monsoon retreat is similarly defined but using the minimum of this sum of wavelet transform coefficients. The use of this method is illustrated using multiple observational datasets and climate model output. The method is compared with existing methods to find onset and retreat dates in three monsoon regions. The method has a similar performance to existing methods that use precipitation thresholds in the North American Monsoon, as shown by the anomalies of precipitation, wind and geopotential anomalies around the onset and retreat dates. The spatial distribution of monsoon onset and retreat in this region was found to be sensibly captured by the wavelet algorithm, illustrating the earlier onset in central western Mexico and the later onset in northwestern Mexico, Arizona and New Mexico. The spatial distribution was shown to be consistent amongst the TRMM, CHIRPS and reanalysis data, suggesting that the method produces similar results in datasets with different climatologies even at the grid-box scale.
The WT method also compares well to a hydrologically defined index (HOWI) in the Indian Monsoon, although the WT better captures the precipitation variations and HOWI the moisture transport. However, the WT method is also able to capture strong differences in moisture transport around the onset and retreat dates, in both models and observations. The WT method obtains a later onset and earlier retreat than in the HOWI index, possibly associated with a lag between the moisture transport in the Arabian Sea (as diagnosed by HOWI) and the actual precipitation over mainland India (as measured by the WT method). The spatial F I G U R E 1 3 Spatial distribution of mean pentads of (a-c) onset, (k-m) retreat as well as (d-f) start and (h-j) end of the MSD in southern Mexico and northern Central America. (a, c, e, g) shows the results for the CHIRPS, (b, d, f, h) for TRMM and (c, f, j, m) for ERA-5 [Colour figure can be viewed at wileyonlinelibrary.com] distribution of onset and retreat dates in the Indian Monsoon region using the WT method seem to be relatively consistent and coherent amongst the observational datasets, as the mean onset date in mainland India was found at pentad 32. Onset is earliest on the western coast of India and the onset date appears to be very homogeneous in central India.
The WT method was extended to characterize the timings and strength of the South American Monsoon bi-modal regime of precipitation and the intervening MSD period, using the same principle but additionally computing the transform over smaller dilations around the onset and retreat dates. By using randomly generated time-series, the spatial distribution of grid-points with a robust MSD signal was found in Cuba, the northwestern coast of Central America and several regions of south and north-eastern Mexico. The MSD in southern Mexico and northern Central America is found to start around pentads 35 and 36 (last week of June) and end around pentad 48 (mid-August) in most observational datasets and reanalysis. To our knowledge, this extension of the WT method provides one of very few methods to characterize the MSD on sub-monthly scales. This method may be potentially useful when diagnosing changes to the characteristics of the MSD in models or observations. Current methods that diagnose monsoon onset and retreat using pentad or daily-mean precipitation timeseries are typically rigid threshold methods. These threshold methods depend on a number of parameters that need to be tuned for a specific monsoon region and for a specific dataset. For instance, the method by G13 used a threshold value specific for the North American Monsoon and specific for the TRMM dataset but also to the limits of the area used for area-averaging the precipitation. In other words, the persistence and threshold values of most of the threshold methods require normalization, statistical treatment or additional tuning to the parameters to account for climatological differences in the datasets which introduces uncertainty. The method by A12 then uses a climatological mean value as the threshold, but in a climate model with a significantly positive bias in the dry winter season of a monsoon this method would be prone to error as the biased seasonal cycle may impose a biased calculation of monsoon onset and retreat.
The WT method is in many ways similar to the agronomical and threshold methods (e.g., Liebmann and Marengo, 2001;Moron and Robertson, 2014), as the implementation of the method uses a subjective determination of the dilation scales; these scales are comparable to the persistence and window parameters of the threshold methods. However, the WT method presented has three main advantages over most threshold methods.
First, the method produces robust results for the Indian and North American monsoon of onset and retreat, and spatial distributions comparable to previous methods (Moron and Robertson, 2014) while not being subject to 'false-hits' nor years without an identification of the onset and retreat dates. In other words, the method provides robust results without requiring further treatment of years with false-hits or undetermined years.
The second advantage of the method is portability, or utility, as the method shows robust and consistent results for three observational datasets, a reanalysis and climate model experiments with varying climate forcing but without any constraint or treatment of the data beforehand and in three different regions with different seasonal cycles. In other words, this method is robust across datasets and regions. In contrast to rigid threshold techniques (e.g., Liebmann and Marengo, 2001), the identification of onset and retreat for each time-series, for example, at each grid-point, is based upon coherent temporal changes within each precipitation time-series while not using parameters determined "a priori" specifically for a region. The WT method then, can be used in any time-series, regardless of the origin of the time-series, without any further change or consideration than those established by the dilations scales determined in Section 2.2.1. The portability of the method also means that the method can be implemented as a "local-scale" method applied at the grid-box scale for high-resolution datasets such as CHIRPS as well as for regional scales using area-averaged time-series.
Third, and in contrast to typical threshold methods, the wavelet method can be applied to climate model output straightforward using the same configuration of dilation scales, a feature of the method that is illustrated by our analysis of several experiments using the Hadley Centre models. The treatment of the data does not require any normalization or statistical treatment even when used for grid-point time-series for different regions or experiments with varying forcing where the seasonal cycle or total annual rainfall may change notably within the model time-series.
While we do not suggest the WT method is better or superior than any others for all purposes, we argue that the method may be potentially useful for multiple purposes due to its portability and robustness. For instance, for climate model inter-comparison analyses which require the determination of onset and retreat in data that may have different seasonal cycles, magnitudes and resolutions. In addition to CMIP-style assessments, the method may be used in analyses of the observed trends in mean onset and retreat dates which may benefit from an independent method of validating the results from other methods or may use several datasets to support their findings. In conclusion, the WT method is designed to be multi-purpose, to be used in varying contexts, from a local grid-scale analysis of the 'agronomical' onset or in global-scale analyses of the onset and retreat of the global monsoon while also being useful for the determination of bimodal regimes, such as those found in Central America and the Caribbean.