Cyclic Interannual Variation in Monsoon Onset and Rainfall in South-Central Arizona, U.S.A.

1. Introduction

1.1. Saguaro Phenology

The association of plant phenological response to the annual cycle of the North American Monsoon (NAM) has been the subject of much study by plant ecologists [1,2,3], but many details are not well characterized. Numerous plants of the Sonoran Desert of North America are partly or wholly dependent on NAM rainfall, but when does monsoon rainfall arrive at any one location? How does arrival timing vary over time and across plant generations? Or does it vary? Appearing as it does in the otherwise extremely hot summer months of June-September, the monsoon rains are vital to plant growth and reproduction [4,5,6]. One plant uniquely dependent on the monsoon is the saguaro cactus (Carnegiea gigantea) which flowers in April, produces fruits from May-June, and disperses seeds by mid-July [7,8,9]. Timing of seed dispersal anticipates the arrival of the first NAM rainfall events in July. Seed germination and seedling establishment requires sufficient NAM rainfalls to trigger germination and then provide moisture, elevated humidity, and lower temperatures to ensure seedling establishment. Seeds apparently do not survive for more than a few weeks after dispersal [10]. Many other species produce and disperse propagules immediately prior to NAM onset, however, as far as is known, at least a portion of those seeds survive the summer and persist in the soil for at least one more season [11,12]. Saguaro synchronization of phenological events with rainfall events and extended periods of favorable climatic conditions is crucial for reproduction. Phenological synchrony with climatic cues is an area of intense ecological research [13,14,15,16].

1.2. North American Monsoon

The synoptic NAM rainfall system emerges from seasonal variations in latitudinal migrations of the Intertropical Convergence Zone (ITCZ, [17,18]). As the ITCZ migrates northward from the equator toward Central America in the northern hemisphere spring, winds along the Pacific coast of Central America and Mexico shift from westerly to southerly transporting increasing quantities of low-level moisture north and northwest along the Mexican coast and eventually entering into Arizona [19,20]. The ITCZ linkage is a relatively recent interpretation made possible primarily through the analysis of global precipitation revealed by satellite sensors [21]. Over the past 30 years, climatologists have achieved considerable success in explaining the causes and mechanisms driving the NAM phenomenon. In particular, the North American Monsoon Experiment (NAME) project of the early 2000s [22]. In addition to compiling impressive datasets from a variety of sensor sources, the NAME researchers worked hard to resolve disagreements regarding the source(s) of NAM moisture and the role of the Gulf of California (GoC) as a mechanistic component of the NAM system [23,24].

Previous research found some indication that the onset and retreat of monsoon weather has changed, perhaps reflecting changes in the ITCZ. Arias, et al. [25] found that the NAM season has shortened due to earlier retreat while the SAM (South American Monsoon) has lengthened by a correlated amount. Ashfaq, et al. [26], using multiple regional projections, predicted broad monsoon onset delays and less overall precipitation, which, if realized, could have significant effects on saguaro phenology.

1.3. Defining Monsoon Onset

The goal of this study was to determine whether saguaro reproductive phenology is synchronized with important climatic events, particularly regarding the timing of seed dispersal and the onset of the monsoon. This first required an objective determination of the date of monsoon onset for a period long enough to discern trends and patterns. We then sought monsoon onset criteria that could be identified from a daily rainfall record, available for many stations in our area from online sources, would yield a single day date, and considered the needs of the plants of the desert. Several methods for identifying monsoon onset were reviewed.

Brenner [27] noted a variety of measurements to indicate that the northward surge of the seasonal maritime airmass from GoC to Arizona has begun, including hourly observations of temperature, dew point, and pressure. In the 1980s meteorologists associated with the U.S. National Weather Forecasting Offices in Phoenix and Tucson began to rely on dew point temperature criteria [28]. Bombardi, et al. [29] catalogued detection methods or criteria from 32 published references, about half of which were solely or partly based on precipitation criteria.

Bowers, et al. [30] found that precipitation amounts triggering seedling emergence in 15 perennial plant species ranged from 17-36 mm. This threshold is higher than necessary to trigger flowering in established perennial plants. Bowers and Dimmitt [2] found that a rainfall event of as little as 9 mm was sufficient to trigger flowering in triangleleaf bursage (Ambrosia deltoidea) and Bowers [8] found that as little as 5 mm of rain triggered flowering in the saguaro cactus. Therefore, for the purposes of this analysis it was determined that the first single-day rainfall of ≥10 mm, from 1 June to 30 September, would be used. This amount is large enough to potentially trigger plant growth and almost every weather station in the dataset recorded on at least one day at some point during the NAM season.

As plant ecologists, we are interested in the influence of environmental signals on phenological triggers and how this might affect the viability of the saguaro cactus and many other plant species. During investigations of the relationship between saguaro phenology and climate, we found interesting patterns in NAM early and late onsets. The purpose of this article is to report the results of our preliminary analyses suggesting previously unreported cyclic interannual variation of in NAM onset and precipitation in south-central Arizona. In subsequent publications we plan to present analyses of the influence of changes in the date of onset of NAM on saguaro phenology.

2. Methods

Precipitation Dataset

The monsoon season was considered to include the period from 1 June to 30 September and monsoon onset was defined as the first day from 1 June to 30 September with single-day rainfall of ≥10 mm. Our primary dataset was constructed with daily precipitation gathered from the Pima County, Arizona, ALERT network of remote reporting rain gauges [31], the Global Historical Climatology Network (GHCN-Daily [32,33]) and Remote Automated Weather Stations (RAWS [34]), most with records from 1990-2022 (Table 1, Figure 1). Data for the full 33-year period from 1990-2022 were obtained from 22 of 33 weather stations. Years of record for the remaining 11 stations ranged from 22 to 32 years.

Variables of interest, which will be referred to as below for the remainder of this paper, extracted from the precipitation dataset included:

• Day-of-year (DOY onset) of the first one-day rainfall events of ≥10 mm from 1 June -30 Sept., 1990-2022. To mitigate possible calendar bias in our method of determining the day-of year of monsoon onset the dates of monsoon onset were counted from the date of the vernal equinox of that year instead of from the first of the year. The date of the vernal equinoxes was obtained from the seasons calculator at timeanddate.com for Tucson, Arizona [35].
• Total monsoon rainfall (MR, mm, 1 June- 30 Sept.).
• Total annual rainfall (AR, mm).
• Proportion of mean annual station rainfall falling from 1 June – 30 Sept (PAR). Each year value in this vector is calculated as the difference between monsoon rainfall of the year and station mean annual rainfall across years.

Of the 33 stations used in the analysis all but three were within a 65 km radius of downtown Tucson, representing 2278 m of elevation from 512 m, at the Organ Pipe Cactus National Monument Headquarters, to 2790 m at the summit of Mt. Lemmon, Santa Catalina Mts. The geographical area covered by the selected stations, 17225 km², was chosen mainly because of the predominance of the Pima County ALERT network in the dataset within the range of the saguaro, but with the addition of stations with sufficient periods of record to sample the region west and south of Tucson, and a few RAWS stations to fill in elevation gaps. The Spatial Autocorrelation tool in ArcGIS Pro v. 3.2 [36] was used to determine whether the 33-year station means of the day after the vernal equinox of the first day rainfall ≥10 mm was influenced by proximity. The input data were weighted to include station elevation (z coordinates) calculated using inverse Euclidean distances.

When required for particular analyses, missing values were filled using the imputeTS package moving average algorithm (R package, ImputeTS [37]). On the few occasions with no rainfall ≥10 mm was recorded during that period (17 of 1,049 DOY records) the highest one-day amount during the period was substituted. The substitution procedure may have led to a slight bias to later DOY onset as at some of these stations the highest amount recorded from 1 June-30 Sept. occurred very late (as late as 30 September).

Least-squares linear regression models (LRMs) were created for the four variables time series listed above as well as (1) DOY onset and MR on elevation, (3-5) MR, AR, and PAR with DOY onset as predictor using base R ([38]). To investigate correlations between stations, Spearman’s Rank statistics for the 33 stations with respect to DOY onset were calculated, also in base R ([38]). Autocorrelation function plots and residuals plots were created for the DOY onset time series using the R forecast package [39]. Normality was investigated using probability density and QQ plots and Shapiro-Wilk tests of the per-station onset time series base R [38]. The time series of the four key variables listed above were tested using multivariate (multi-site) Mann-Kendall monotonic trend (R package trend, [40]). Cross-correlation tests were run with MR, AR, and PAR time series as response and DOY onset as the predictor using base R [38].

Patterns in scatterplots for each variable (DOY onset, MR, AR, PAR) suggested non-linearity that generalized additive models (GAMs) might elucidate. GAMs are widely employed in many fields when non-parametric analyses and relaxation of normality assumptions are desired [41,42,43]. Fit curves with basis dimensions set to k = 10 were generated from GAMs developed in R using the MGCV package [44] for the four variables. The results were useful in estimating sinusoidal wave periods and coefficient starting points for sinusoidal regression models (SRMs) using the MATLAB Curve Fitter tool [45] of the form:

f(y)=a sin(bx+c)+d

where y = DOY onset, MR, AR, or PAR, a = amplitude, x= year, b = change in period, c = horizontal shift, and d = vertical shift (calculated in base R [38].

3. Results & Discussion

3.1. Effect of Elevation on DOY Onset & Precipitation

In our NAM study area, orography strongly influences precipitation and at most sites there is a significant positive relationship between monsoon precipitation and elevation [46,47]. Figure 2 illustrates this with LRM of mean station MR on station elevation (R²=0.891, p-value=<.001). In addition, DOY onset data regressed on elevation (Figure 3) also suggested a positive relationship (R²=0.540, p-value=<.001). Thus, an early onset is associated with higher MR along the elevational gradient. Higher elevation sites generally have a greater probability of rainfall ≥10 mm on any day (both from the increased frequency of rainfall events and from the greater precipitation received from each event) and are, therefore, more likely to achieve our rainfall threshold (≥10 mm/day) sooner rather than later in any arbitrarily chosen period, including the monsoon season. The association of earlier onset with increasing elevation in our region is not unexpected, however, it is surprising that our DOY onset data, mixing untransformed gauge data from a broad elevational range should nevertheless exhibit the coherent cyclic pattern of interannual variation described in the next section. Early in our analysis we removed the effects of elevation on DOY onset variance, however, this had almost no effect on subsequent analyses and was not pursued. The positive correlation between monsoon precipitation and elevation observed in our gauge ensemble was not observed in other areas, including some relatively close by, such as northwestern Mexico, where distance from the Pacific Ocean coastline and the Sierra Madre Occidental are more reliable predictors ([48,49].

3.2. Monsoon Day-of-Year Onset

Table 2 lists the 33 weather stations, mean DOY onset, MR, AR, and PAR for each station. (Note that DOY onset is counted from the day of vernal equinox of each year). The mean DOY onset of all 33 stations from 1990-2022 was 16 July (min.=20 June, max.=8 August), 117.6 days from the vernal equinox date (usually March 20). Figure 4 is a correlation matrix (Spearman’s Rank) of DOY onset for the 33 weather stations listed in Table 2 sorted by elevation. Of the 528 correlations, 346 are significant (p-value ≤0.05), all are positive correlations, and all but seven of those are moderately correlated (ρ < 0.75). The mean value of the 528 correlation coefficients is 0.40. The ArcGIS Pro v. 3.2 Spatial Autocorrelation tool [36] of DOY onset means for each station produced Moran’s I = -0.005, z = 0.318, p-value = 0.750, indicating not different from random. ArcGIS Pro v. 3.2 Spatial Nearest Neighbor tool [36] indicates spatial distribution not different from random (p-value=0.654) and the mean between-station distance is 13.5 km. At nearest neighbor distances of 17.0 km and 11.3 km Gebremichael, et al. [50], Mascaro, et al. [51], respectively, also found low between-station correlations.

Indeed, although DOY onset is associated with elevation, visual inspection of Figure 4 suggests that sites at similar elevations do not appear more likely to be correlated with respect to each other. Even very close stations may exhibit low correlations. Sabino Dam and Saguaro (Table 1), separated by 337 m horizontal and 56 m elevation, are correlated with ρ = 0.385, while Organ Pipe Cactus National Monument and Mt. Lemmon, ρ = 0.431, are separated by 197 km and 2274 m elevation.

The partial independence of monsoon rainfall gauge readings in the desert and semidesert lands of southwestern North America has long been recognized ([52]), which is how our gauge-based network is treated for the purpose of this analysis.

Figure 5 depicts the probability density plot and normal curve of the DOY onset dataset (all 33 stations, 1990-2022, 1,049 data points) showing a left-skew reflecting the fewer, but more tightly clustered early onset DOY onsets characteristic of the sinusoidal wave troughs. Shapiro-Wilk tests of the per-station onset time series indicate non-normality in the onset day for two-thirds of stations and non-normality of onset distribution when the data are tested per-year in every year from 1990-2022. Monsoon DOY onset is therefore likely to be normally distributed for the station when accumulated over time, but non-normal for DOY onset for any given year across the stations.

We suggest that the skewed DOY onset timing is the result of at least two conditions: (1) an imbalanced effect of the fewer early DOY onset years that exhibit less variance as discussed in the next section, and (2) the beginning date of the monsoon was set at 1 June while the end date was set at 30 September, a span of 122 days. Mean DOY onset date was 16 July, 46 days from 1 June and 76 days from 30 September. There is more time after the mean DOY onset date for stations to achieve the ≥10 mm criterion than before.

3.3. GAMs, LRMs, and SRMs

Figure 6, Figure 7, Figure 8 and Figure 9 present scatterplots, GAMs, LRMs, and SRMs for DOY onset, MR, AR, and PAR (1 June-30 September) for the 33 weather stations listed in Table 2. Summary statistics are listed in Table 3 and crest-trough periods and amplitudes derived from the GAM curves are listed in Table 4.

Notwithstanding the low p-values listed in Table 3 for the linear regressions of DOY onset, MR, AR, and PAR shown in Figure 6, Figure 7, Figure 8 and Figure 9, low R² values indicate very little variance, if any, is explained. This was due, at least in part, to the confounding effects of the sinusoidal wave pattern. Tests of the time series with multivariate (multi-site) Mann-Kendall monotonic trend (R package, trend, [40]) showed no trend. Univariate Mann-Kendall monotonic trend (R package, trend, [40]) were run on the individual station time series finding no trend in 28 of the 33 stations and of the five with p-value ≤0.05, all had negative trends.

The four plots of GAMs presented in Figure 6, Figure 7, Figure 8 and Figure 9 and Table 3 were run with basis dimensions, k = 10. Increasing basis dimensions (up to 33) may explain as much as 50% of deviance but, in our opinion, over-explains the data. At k=10, GAMs of DOY onset and PAR explain 19.1% and 21.4% of deviance, respectively, the MR GAM, explains 13.3%, and the AR GAM explains 5.4%.

The three GAMs and SRMs shown in Figure 7, Figure 8 and Figure 9 are more or less faithful representations of the GAM and SRM fit curves in Figure 6 but in mirror image. Nevertheless, time series autocorrelation analyses (R package, trend, [40]) of DOY onset failed to find correlations with AR at any lag other than +2. However, linear models using the means of DOY as predictor and MR and PAR as responsors produced R² of 0.431 and 0.418, respectively (p-value = <0.001). As have many other researchers, we found that early DOY onset produces more rain than late onsets [26,53,54].

The diagnostic bar charts, residual scatterplots, and QQ plots provided in Figure 6, Figure 7, Figure 8 and Figure 9, all show the same left-skew, and slight to substantial departures from model expectation. Slight trends in the LRM lines for DOY onset, MR, and PAR, while presenting below threshold (<0.05) p-values, explain very little variance. Taken together however, they suggest the possibility that DOY onset is advancing in the year, MR is increasing, and PAR is increasing, implying a slight decrease in October-May rainfall.

Scatterplots of the time series of DOY onset, MR, AR, and PAR received during the monsoon suggest that, with respect to DOY onset, the sinusoidal troughs exhibit less variance than the crests. The uneven distribution of variance between the troughs and crests discussed in the previous section may be partly responsible for the low deviance and R² values produced in our GAM analysis and SRM fit.

The sinusoidal wave pattern seen in the scatterplots is reversed in the GAMs and SRMs of monsoon precipitation, annual precipitation, and the proportion of annual precipitation falling during the monsoon. Similar to the DOY onset waves, but in mirror image, these waves are broad in the troughs and narrow in the crests reflecting the predominance of late onsets and less rainfall. The unbalanced distribution of variance coupled with the unequal spans of crests and troughs may explain the fit statistics observed in the GAMs and SRMs.

It is noteworthy that the NAME project previously referenced focused on the monsoon season of I June – 30 September 2004 [22], coinciding with an onset crest (Figure 7a) and monsoon precipitation trough (Figure 8a).

Table 4 is presented recognizing that the specific shapes of the curves rendered in the GAMs are partly subject to the choice of basis dimensions (k=10). Nevertheless, although period and amplitude estimates of the GAM curves (Figure 6, Figure 7, Figure 8 and Figure 9) were used as starting points, each SRM required multiple iterations to arrive at the model selected and these match the GAMs relatively well. The DOY onset curve period, 3140 days (8.6 years) is calculated from the four distinct crests in the GAM curve (Figure 6), while the MR and AR and the PAR curves include only three crests in the 1990-2022 dataset, 3039, 3050, and 2969 days (8.3, 8.4, and 8.1 years) respectively. Note that the DOY onset wave period and the periods of the MR, AR, and PAR are not identical. The waves may be out of synchrony or, more likely, the data may not be sufficiently robust.

The sinusoid wave amplitudes (crest minus trough) produced by the GAMs are very substantial. The smallest amplitude in the DOY onset time series GAM is 14 days, the largest, 29 days or more than four weeks. The mean amplitudes of MR, AR, and PAR are 90 mm, 42 mm, and 25% respectively.

3.4. Previous Research

Climatologists have been researching the nature of the NAM for more than a century [55]. In assessing our results, we surveyed the relevant literature, a portion of which is listed in Table 5. The aims, priorities, and tools available to researchers over time have evolved from recognition, to characterization, investigations of origin, modelling, linkage to global systems, and, most recently, the future of the monsoon in a warming climate. These focus areas overlap and continue. Throughout this extensive history of research the identification and analysis of monsoon onset and retreat has relied on several approaches: single day rainfall thresholds, wind measurements, multi-day dewpoint and rainfall thresholds, combinations of climate variables, rainfall indices, and more.

Table 5 lists 21 NAM research contributions, only two of which, Ellis, et al. [56], García-Franco, et al. [57], were devoted to identifying dates of monsoon onset. Of the 20, two authors used multiple onset identifiers [27,58], 11 relied on variously constrained multi-day precipitation thresholds [25,26,54,59,60,61,62,63,64,65,66,67], five used abrupt or anomalous changes in rainfall amounts [27,53,58,68,69], three favored a multi-day dewpoint threshold [27,56,70], three prioritized change in wind direction [27,55,58], one employed multiple surface measurements [27], and one developed a novel wavelet transform procedure [57]. Each author chose a monsoon onset definition that would facilitate analysis or support a research goal as did we in choosing our own unique definition. Our interest lay in the need for desert plants at the tail end of a lengthy foresummer drought amid increasingly harsh temperatures to take advantage of the coming rains as quickly as possible. The ≥10 mm criterion was chosen as the least amount that would potentially trigger plant growth to support investigations of the phenology of desert plants, not monsoon season climatology. As it happened, our approach may have revealed an important pattern in DOY onset and precipitation that has so far escaped notice.

With one exception ([48]) among all research studies cited in this paper, analyses were performed using summarized data, usually by extracting and analyzing mean values, often in gridded datasets or compiled as means of variables from terrestrial stations. This was presumably done to avoid violations of independence of observations. Our dataset was extracted from terrestrial sensors and the main analysis was performed with daily rainfall data (although we report mean values) from which we extracted single-day yearly onsets for each station for each year 1990-2022. As reported above, we tested the assumption of independence in our DOY onset variable station means using the ESRI ArcGIS spatial autocorrelation tool [36] and found the resulting Moran’s I indicated not different from random. The sinusoidal wave pattern observed in the scatterplot is likely to be real and worthy of further investigation.

3.5. Saguaro Phenology

The research described in the previous section suggests sinusoidal waves in timing of DOY onset, MR, AR, and PAR from 1990-2022. It is unknown whether that pattern extends back in time to any extent. Saguaros are very long-lived, up to 250 years [71], and it would be very interesting to know whether, or to what degree saguaros have adapted to regular cycles of abundance and shortage. Even if the sinusoidal wave only appeared de novo in 1990, saguaro reproduction, which as we have noted, is highly dependent on the monsoon as it is now. As noted above, the variations in DOY onset and MR reported here are quite substantial. The sinusoidal waves have amplitudes of up to 30 days and nearly 100 mm, respectively.

The possible effect of changes in DOY onset timing and monsoon precipitation, over time on saguaro phenology bears further research in light of predicted changes in saguaro distribution under the influence of anthropogenic climate change. Albuquerque, et al. [72] predict loss of currently occupied saguaro habitat in the next half century and Renzi, et al. [7], using saguaro phenological data from a population of 151 saguaros monitored from 2004-2013 found that saguaro flowering was increasingly delayed while the number of flowers produced declined overall. These reports do not, at least on the surface, appear to reflect any relationship with the sinusoidal waves suggested in our data. Attempts have been made to analyze multi-decadal patterns of recruitment and growth [73], but it is highly unlikely that the very long-lived saguaro cactus would adapt to this short-period cyclic climatology.

Table 5. List of monsoon onset papers, datasets, and approaches. Results of selected literature survey of synoptic NAM characteristics addressing onset, duration, or precipitation amounts, solely or in combination with other climatological investigations.

	Reference	Aim of study	Onset identifier(s)	Geography
1914	Huntington [55]	Describe monsoon climate of Ariz., New Mex.	Change in wind direction from westerly to southerly.	Ariz., New Mex.
1955	Bryson and Lowry [68]	Map shift in dominant air masses signaling monsoon season.	Sharp changes in Raininess Index.	Ariz.
1973	Brenner [27]	Describe antecedent monsoon climate of Ariz., investigate possible GoC connection.	Changes in several surface observations.	Ariz.
1993	Douglas, et al. [58]	Describe synoptic climatology Mex., U.S. monsoon, identify moisture source(s).	Rainfall amounts and wind shifts.	Mex., SW U.S.
1997	Higgins, et al. [59]	Diagnose atmospheric conditions preceding monsoon onset.	Precip. index magnitude & duration, ⫹0.5 mm day⫺1 for 3 days.	Ariz., western New Mex.
1998	Higgins, et al. [60]	Extend earlier work, diagnose variability of Mex.-U.S. summer precip.	Precip. index magnitude & duration, ⫹0.5 mm day⫺1 for 3 days.	Ariz., western New Mex.
1999	Higgins, et al. [61]	Extend earlier work, identify factors influencing variability.	Varies by region; precip. index magnitude & duration, ⫹0.5-2 mm day⫺1 for 3-5 days.	Ariz., western New Mex., SW Mex.
2002	Mitchell, et al. [62]	Onset & evolution of monsoon rainfall related to GOC SSTs.	Varies by region; precip. index magnitude & duration, ⫹0.5-2 mm day⫺1 for 3-5 days..	GoC, Ariz., New Mex.
2004	Ellis, et al. [56]	Develop regional criteria for monsoon onset.	Mean daily dew-point threshold 12.2º or 12.8º for 3 days.	SW U.S.
2007	Grantz, et al. [69]	NAM timing, rainfall amount, large-scale trends drivers.	Percentile of monsoon precip. threshold.	Ariz., New Mex.
2008	Liebmann, et al. [53]	Monsoon onset dates Mex. & U.S. climatology variability, season length, rate, SST topography influences.	Anomalous rainfall accumulation.	Mex., SW U.S.
2009	Turrent and Cavazos [63]	Clarify monsoon forcing with land-sea thermal flux and influence on onset.	First 5-day period mean core region precip. >1 mm.	NW Mex.
2011	Crimmins, et al. [70]	Relate summer flowering onset to monsoon climatological events.	3-day mean daily dewpoint threshold at TUS.	Finger Rock Canyon, Pima Co., Ariz.
2012	Arias, et al. [64]	Determine multidecadal variations in monsoon seasonality & strength.	Pentad after which the rain rate > than annual mean rain rate in 6 of 8 preceding pentads and after which rain rate > than annual mean in 6 of 8 pentads.	NW Mex.
2015	Arias, et al. [25]	Investigate climate dynamics influencing variability of NAM and SAM.	Pentad after which the rain rate > than annual mean rain rate in 6 of 8 preceding pentads and after which rain rate > than annual mean in 6 of 8 pentads.	North America, South America
2017	Meyer and Jin [54]	Correct regional and global models for projected future climate effects on monsoon dynamics.	First occurrence after 1 May with three consecutive days with at least 0.5 mm day−1.	SW U.S., Mex.
2021	García-Franco, et al. [57]	Introduce new onset identifier.	Wavelet transform (maximum sum of wavelet coefficients) on precip. time series.	Mex., India
2021	Fonseca-Hernandez, et al. [65]	Analyze mixing mechanisms responsible for temporal and spatial variations of GoC boundary layer during NAM onset.	First day first sequence of five consecutive days mean precip. rate =/> 2 mm/day.	NW Mex.
2021	Ashfaq, et al. [26]	Construct regional climate model of monsoon change with increased greenhouse gas forcing.	Pentad after minimum seasonality of precip. value, similar to Bombardi and Carvalho [74].	Global monsoon domains
2024	Duan, et al. [67]	Identify and analyze NAM extreme events.	First 1 mm/day for 5 days after 1 June.	GoC and surrounding lands

Table 5. List of monsoon onset papers, datasets, and approaches. Results of selected literature survey of synoptic NAM characteristics addressing onset, duration, or precipitation amounts, solely or in combination with other climatological investigations.

	Reference	Aim of study	Onset identifier(s)	Geography
1914	Huntington [55]	Describe monsoon climate of Ariz., New Mex.	Change in wind direction from westerly to southerly.	Ariz., New Mex.
1955	Bryson and Lowry [68]	Map shift in dominant air masses signaling monsoon season.	Sharp changes in Raininess Index.	Ariz.
1973	Brenner [27]	Describe antecedent monsoon climate of Ariz., investigate possible GoC connection.	Changes in several surface observations.	Ariz.
1993	Douglas, et al. [58]	Describe synoptic climatology Mex., U.S. monsoon, identify moisture source(s).	Rainfall amounts and wind shifts.	Mex., SW U.S.
1997	Higgins, et al. [59]	Diagnose atmospheric conditions preceding monsoon onset.	Precip. index magnitude & duration, ⫹0.5 mm day⫺1 for 3 days.	Ariz., western New Mex.
1998	Higgins, et al. [60]	Extend earlier work, diagnose variability of Mex.-U.S. summer precip.	Precip. index magnitude & duration, ⫹0.5 mm day⫺1 for 3 days.	Ariz., western New Mex.
1999	Higgins, et al. [61]	Extend earlier work, identify factors influencing variability.	Varies by region; precip. index magnitude & duration, ⫹0.5-2 mm day⫺1 for 3-5 days.	Ariz., western New Mex., SW Mex.
2002	Mitchell, et al. [62]	Onset & evolution of monsoon rainfall related to GOC SSTs.	Varies by region; precip. index magnitude & duration, ⫹0.5-2 mm day⫺1 for 3-5 days..	GoC, Ariz., New Mex.
2004	Ellis, et al. [56]	Develop regional criteria for monsoon onset.	Mean daily dew-point threshold 12.2º or 12.8º for 3 days.	SW U.S.
2007	Grantz, et al. [69]	NAM timing, rainfall amount, large-scale trends drivers.	Percentile of monsoon precip. threshold.	Ariz., New Mex.
2008	Liebmann, et al. [53]	Monsoon onset dates Mex. & U.S. climatology variability, season length, rate, SST topography influences.	Anomalous rainfall accumulation.	Mex., SW U.S.
2009	Turrent and Cavazos [63]	Clarify monsoon forcing with land-sea thermal flux and influence on onset.	First 5-day period mean core region precip. >1 mm.	NW Mex.
2011	Crimmins, et al. [70]	Relate summer flowering onset to monsoon climatological events.	3-day mean daily dewpoint threshold at TUS.	Finger Rock Canyon, Pima Co., Ariz.
2012	Arias, et al. [64]	Determine multidecadal variations in monsoon seasonality & strength.	Pentad after which the rain rate > than annual mean rain rate in 6 of 8 preceding pentads and after which rain rate > than annual mean in 6 of 8 pentads.	NW Mex.
2015	Arias, et al. [25]	Investigate climate dynamics influencing variability of NAM and SAM.	Pentad after which the rain rate > than annual mean rain rate in 6 of 8 preceding pentads and after which rain rate > than annual mean in 6 of 8 pentads.	North America, South America
2017	Meyer and Jin [54]	Correct regional and global models for projected future climate effects on monsoon dynamics.	First occurrence after 1 May with three consecutive days with at least 0.5 mm day−1.	SW U.S., Mex.
2021	García-Franco, et al. [57]	Introduce new onset identifier.	Wavelet transform (maximum sum of wavelet coefficients) on precip. time series.	Mex., India
2021	Fonseca-Hernandez, et al. [65]	Analyze mixing mechanisms responsible for temporal and spatial variations of GoC boundary layer during NAM onset.	First day first sequence of five consecutive days mean precip. rate =/> 2 mm/day.	NW Mex.
2021	Ashfaq, et al. [26]	Construct regional climate model of monsoon change with increased greenhouse gas forcing.	Pentad after minimum seasonality of precip. value, similar to Bombardi and Carvalho [74].	Global monsoon domains
2024	Duan, et al. [67]	Identify and analyze NAM extreme events.	First 1 mm/day for 5 days after 1 June.	GoC and surrounding lands

4. Conclusions

The dependency of saguaro population biology on monsoon precipitation prompted us to gather precipitation data from 33 weather stations in south-central Arizona and identify the date of monsoon onset at each station from 1990-2022. In the near future we plan to update the saguaro phenology dataset in Renzi, et al. [7] to address this question.

The results reported here suggest that the DOY onset time series exhibits a sinusoidal wave fit to a GAM with period of 8.6 years and mean crest-trough amplitude of 20 days. The relationship is weak, explaining only 19.1% of deviance. We suggest some preliminary explanations for greater data variance from troughs to crests.

We report interesting and, we think, very important patterns in selected NAM characteristics. We can offer no plausible explanation except to suggest interplay between sub-global climate systems that leads to spatially recurring systems. Before the drivers of the monsoon patterns suggested here can be elucidated much work remains to be done to establish: (1) whether, or which, of the patterns, especially timing of DOY onset, are real, (2) what is the spatial extent of the patterns, (3) how long before 1990 do the patterns extend. and (4) can we create models sufficiently robust to support forecasting?

It is hoped that future work on the DOY onset features found in our data will be explained and clarified with greater resolution than is presented here. We look forward to insights into the nature of the sinusoidal wave in the DOY onset data, clarification of the geospatial characteristics of DOY onset timing, and we hope that some progress may be made in extending analysis of DOY onset to substantially before 1990.