Introduction

In recent decades, one of the most damaging agroclimatic phenomena worldwide is the increase in the trend of meteorological droughts (MeD) [1-3]. Clarke [4], Spinoni et al. [5], and Syed et al. [6] state that the significant trend of MeD can be studied using linear regression (SLT) and the non-parametric tests (SNT) of Mann–Kendall [3, 7-14]. These MeD trends have been applied to the agricultural standardized precipitation index (aSPI) [6] and reconnaissance drought index (RDI) [9, 10]), mainly due to the fact that they only require effective precipitation (fP) as input data [15-17] and total precipitation and potential evapotranspiration (rP–PET) [16, 18, 19]. Specifically, Syed et al. [6] and Merabti et al. [10] showed that the higher temporal variability of aSPI and RDI in Saudi Arabia and Algeria respectively can be explained by the first observed principal component (Z PC–1o). Additionally, Syed et al. [6] observed that the trend of aSPI and RDI in Saudi Arabia is highly associated with the variation of the Pacific decadal oscillation (PDO). On the other hand, Zhang et al. [20] and Li et al. [ 21] found in China that the PDO and Atlantic multidecadal oscillation (AMO) function as predictors of MeD. According to Méndez and Magaña [22] and Llanes et al. [23], the MeD in northern Mexico can be associated with and predicted through the –PDO and +AMO phase anomalies since they promote the occurrence of La Niña events (shortage of rP and PET). The authors of this study have made various contributions in this field of investigation in Sinaloa; for example, Llanes et al. [14] calculated the SNT of MeD using the standardized precipitation evapotranspiration index (the trend of MeD was not modeled). Llanes et al. [15] found that of four indices of MeD that use fP (aSPI and effective reconnaissance drought index) and rP (standardized precipitation index and RDI) as explanatory variables, aSPI and RDI are respectively the indices with the greatest sensitivity to predict rainfed maize yield, which is the crop with the largest sown area in Sinaloa. Finally, Llanes [16] used PDO–AMO to predict aSPI–RDI (without evaluating the trend).

Therefore, this study is a pioneer for Sinaloa in predicting SLT–SNT of aSPI–RDI (dependent variables) using PDO–AMO (independent variables).

The goal is to model the trend of MeD using the PDO– AMO indices. This study contributes to the prediction of MeD trends in one of the states with the highest agricultural production in Mexico [15]. Understanding these trends can help ensure regional and national food sovereignty, especially in abnormally dry periods [22].

Material and Methods

Generalities

In this study (Sinaloa state) and for the period 1981– 2017, the daily data of rP, minimum (minT) and maximum temperatures (maxT) were obtained from Comisión Nacional del Agua–Servicio Meteorológico Nacional (CONAGUA–SMN) [24]. SLT–SNT were determined for 12 weather stations after calculation of aSPI and RDI (using the DrinC program). To obtain the greatest variability of SLT–SNT, Z PC–1o was extracted. The monthly series of PDO and AMO for the period 1981–2017 were obtained from the National Oceanic and Atmospheric Administration (NOAA) [25]. To find the relationship between SLT–SNT and PDO–AMO, a Pearson correlation (Pe) for Z PC–1o vs PDO–AMO was applied.

For the calculation of SLT–SNT and the correlation of Z PC–1o vs PDO–AMO, the following scales and time steps were employed:

  • 1) Nov–Jan, as it is the quarter with the highest sensitivity for predicting rainfed maize yield, considering aSPI and RDI as independent variables [16].

  • 2) Jun–Nov, since according to Llanes et al. [26], it is the semester where the rP with the highest annual contribution occurs.

  • 3) Jul–Sept and Mar–Aug, for being respectively the quarter and semester with the highest annual contribution of rP and PET. In addition, in these two-time steps, PDO–AMO indices have more explanatory power for predicting aSPI and RDI [16].

  • 4) Oct–Sept was considered as the agricultural year, in addition to being the order of data series in DrinC software, to calculate aSPI and RDI [27].

After calculating the correlation of Z PC–1o vs PDO– AMO, in order to model the first calculated principal component (Z PC–1c), a linear regression was applied to Z PC–1c (dependent variable) vs [Z PDO–6 (Mar– Aug), Z AMO–3 (Jul–Sep) and Z AMO–6 (Mar–Aug)] (independent variables) [28]. To validate the predictive sensitivity for the trend of Z PC–1c, a Pe was applied between Z PC–1o and Z PC–1c for the period 1981–2017. After a standardized normalization [29] and a Shapiro– Wilk normality analysis [27, 28] was done, the trend of Z PC–1c was calculated using linear regression). The correlations Z PC–1o vs Z PC–1c and the linear trend of Z PC–1c were significantly different from zero.

Study Area

This study was carried out in the state of Sinaloa in the northwest of Mexico (Fig. 1). This state is known worldwide for its high volumes of agricultural production (“the breadbasket of Mexico”). For example, it is the national leader in the production of grains, legumes, and vegetables [30]; however, due to its extreme temperature [30] and precipitation gradients [26], the state is vulnerable to the occurrence of intense PET and MeD [23]. In the data series of this study, the minimum and maximum average values of rP and PET ranged from rP = 439.68 mm (El Playón) to rP = 1478.75 mm (Potrerillos) [31] and from PET = 1291.17 mm (Potrerillos) to PET = 2029.60 mm (Ixpalino). The average aridity index (A.I. = 0.54, dry sub-humid) ranged from A.I. = 0.25 (semi-arid, El Playón) to A.I. = 1.15 (wet, Potrerillos).

Data

For the period 1981–2017, daily data of rP, minT and maxT were obtained from 70 weather stations located in the state of Sinaloa. The data were obtained from the National Water Commission–National Meteorological Service [24] at https://smn.conagua.gob.mx/es/climatologia/informacion-climatologica/informacion-estadistica-climatologica. The proportion of missing data was less than 10% at only twelve stations; therefore, the series from these stations were selected for this study. Missing data were estimated with the multiple imputation method, using the Markov chain Monte Carlo algorithm (MCMC) [32], since according to Aieb et al. [33] and Johnson [34], it is the best choice for imputation when the missing data <10%.

Fig. 1.

Study area. Location of weather stations.

https://www.pjoes.com/f/fulltexts/187162/PolJEnvironStud-2025-34-2-1515-1527_Fig_1_min.jpg

Mathematical Equations to Calculate the Indices of Meteorological Drought (MED)

Agricultural Standardized Precipitation Index (aSPI)

The index aSPI is a modification of the standardized precipitation index [15, 35, 36]; the difference is that instead of using rP, aSPI only uses fP. In this study, fP was calculated with the USDA method–version CROPWAT [37], which, according to Tigkas et al. [15], has been used successfully and on multiple occasions in semi-arid conditions. Equations 1 and 2 characterize fP:

(1)
fP=rP·(125−0.2rP)125;for rP≤250 mm
(2)
fP=(0.1·rP)+125;for rP>250 mm,

where fP is the effective precipitation and rP is the total precipitation.

To calculate aSPI, an index created by Tigkas et al. [15] and applied by Llanes et al. [16] and Javed et al. [38], the following procedure was applied: first, the fP series was fitted to the gamma distribution (Equation 3):

(3)
g(x)=1βαΓ(α)xα−1e−xβ,for x>0,

where the shape and scale parameters are made up of α and β, respectively. Γ(α) denotes the gamma function and x = fp. Equations 4 and 5 describe the maximum likelihood functions respectively.

(4)
α=14A(1+1+4A3)
(5)
β=x¯α

where n is the number of observations (Equation 6) and where

(6)
A=ln(x¯)−∑ln(x)n.

To obtain a cumulative probability G(x) from fP, Equation 7 was used:

(7)
G(x)=∫0xtα−1e−tβdtΓ(α).

Since fP can contain null values, Equation 8 [39] satisfies the cumulative probability:

(8)
H(x)=q+(1−q)·G(x).

Using Equations 9–11, the cumulative distribution was converted to a normal distribution [40] to obtain the aSPI index:

(9)
aSPI=−(t−c0+c1t+c2t21+d1t+d2t2+d3t3),0<H(x)≤0.5where t=ln1H(x)2
(10)
aSPI=t−c0+c1t+c2t21+d1t+d2t2+d3t3,0.5<H(x)≤1.0where t=ln1[1−H(x)2]
(11)
and c0=2.515517,c1=0.802853,c2=0.010328,d1=1.43278,d2=0.189269,d3=0.001308.

Reconnaissance Drought Index (RDI)

To obtain RDI, rP and PET were first calculated. In this study, PET was calculated using the Hargreaves–Samani method [16, 17, 41–46]. Equation 12 defines PET:

(12)
PET=0.0023·(mT+17.8)·(maxT−minT)0.5·Ra,

where PET is the potential evapotranspiration [mm month–1]; mT, maxT and minT are respectively the average, maximum and minimum temperatures [°C month−1]; and Ra is the extraterrestrial solar radiation [mm month–1]. To calculate RDI, which was created by Tsakiris [18] and applied by Syed et al. [6] and Vangelis et al. [46], Equations 13 and 14 were used:

(13)
αk(i)=∑j=1krPij∑j=1kPETij,i=1(1)·N, where j=1(1)·k
(14)
RDIn(i)=ak(i)ak(i)―−1, where ak(i)― is the arithmetic mean of ak(i),

where (ak) of RDI is calculated for the i−th year on a time basis of k (months), and j is the month of the i−th year. N is the total number of years of the available data.

According to the explanation of Vangelis et al. [46], it was assumed that ak satisfactorily follows a lognormal distribution. Then, using Equation 15, the standardized RDI was calculated:

(15)
RDIst(i)=y(i)−y¯σ^y,

where y(i) is ln(ak(i)), y¯ is its arithmetic mean and σ^y is its standard deviation.

Mathematical Equations to Determine Significant Non-Parametric Trends (SNT)

Mann–Kendall Test

Since not all MeD data series presented normality, the SNT were first calculated. Using Equation 16, the Mann– Kendall test was applied [7, 8].

(16)
S=∑i=1n−1∑j=i+1nsgn(xj−xi),

where n is the number of data points in the series, xj and xi are the data values of items j and i; with the condition j>i; and sgn(xj−xi) is the sign function. The sign function is described by Equation 17:

(17)
sgn(xj−xi)={+1,if xj−xi>00,if xj−xi=0−1,if xj−xi<0.

The variance function is defined by Equation 18:

(18)
Var(S)=n(n+1)(2n+5)−∑i=1mti(ti−1)(2ti+5)18,

where n is the number of data points in the series, m is the number of linked groups and ti is the number of ties in measurement i. Since in this study n>10 , the standard normal test Zs was determined with Equation 19:

(19)
Zs={S−1Var(S),if S>00,if S=0S+1Var(S),if S<0.

For Zs>0, it defines an increasing trend and when Zs<0, a decreasing trend is identified. A SNT is defined when Zs>|1.96| [3].

Sen’s Slope Estimate

This statistical test was applied through the following steps: the slope of a trend in the sample from pairs of data was estimated [12] (Equation 20).

(20)
Qi=xj+xkj−k,for i=1,…,N,

where xj and xk are the data values at times j and k, with the condition j>k. In the event of finding only one datum for each period, Equation 21 was applied:

(21)
N=n(n−1)2,

where n is the number of periods. If there were multiple observations in one or more time periods, Equation 22 was applied:

(22)
N<n(n−1)2,

where is the total number of observations. N values from Qi were classified from highest to lowest, and Sen’s slope estimator was calculated with the median, Equation 23 was used:

(23)
Qmed={Q[N+12],if N is oddQ[n2]+Q[N+22],if N is even.

The sign from Qmed was the same as that of Zs (Equation 19) and the magnitude from Qmed was the steepness of the trend (Equation 23).

Mathematical Equations to Determine Significant Linear Trends (SLT)

Line Fitting by the Least Squares Method

After calculations of the SNT, SLT were determined. To calculate SLT and based on what was described by Devore [47], it is known that the simplest deterministic relationship between two variables and is the linear relationship, which is expressed by Equation 24:

(24)
yi=β0+β1xi.

Using Equation 25, the vertical deviation of a point (xi,yi) of a line is defined as:

(25)
y=b0+b1x.

Therefore, the height of the point (xi,yi) can be expressed by Equation 26.

(26)
yi−(b0+b1xi).

Then the sum of the squared vertical deviations of the points (x1,y1),…,(xn,yn), is defined by Equation 27:

(27)
f(b0,b1)=∑i=1n[yi−(b0+b1xi)]2,

where the estimated points β0 and β1, denoted by β^0 and β^1, are called the least squares estimates; that is, they are the values that minimize f(b0,b1). β^0 and β^1 are expressed by Equation 28:

(28)
f(β^0,β^1)≤f(b0,b1).

Equation 28 can be applied to any value of b0 and b1. Therefore, the estimated regression line (least squares line) is denoted by Equation 29:

(29)
y=β^0+β^1x.

The minimizing values of and are found by taking the partial derivatives of f(b0,b1) with respect to b0 and b1 (Eq 28). If the function of a single variable is set equal to zero, the result is Equation 30:

(30)
f′(b)=0.

Equations 31 and 32 are then solved:

(31)
∂f(b0,b1)∂b0=∑i=1n2(yi−b0−b1xi)(−1)=0
(32)
∂f(b0,b1)∂b1=∑2(yi−b0−b1xi)(−xi)=0.

The common factor –2 is canceled and the equations are rearranged, resulting in the system of normal Equations 33 and 34:

(33)
nb0+(∑i=1nxi)b1=∑i=1nyi
(34)
(∑i=1nxi)b0+(∑i=1nxi2)b1=∑i=1nxiyi.

These two equations are linear in the two unknowns b0 and b1. As long as all the values of xi are not identical, the least squares estimate is the only solution of the system. The least squares estimate of the slope coefficient b1 of the true regression line is denoted by Equation 35:

(35)
b1=β^1=∑i=1n(xi−x¯)(yi−y¯)∑i=1n(xi−x¯)2=SxySxx.

The formulas for calculating the numerator and denominator of β^1 result in Equations 36 and 37:

(36)
Sxy=∑i=1nxiyi−(∑i=1nxi)(∑i=1nyi)n
(37)
Sxx=∑i=1nxi2−(∑i=1nxi)2n.

Therefore, the intersection b0 of the true regression line, which was estimated by least squares, is expressed by Equation 38:

(38)
b0=β^0=∑yi−β^1∑xin=y¯−β^1x¯.

To determine SLT, in this study the same methodology was used as in Llanes et al. [26], where the sign and magnitude of the slope coefficient β^1 were considered as the phase and magnitude of the SLT increase, respectively.

Meteorological Drought (Med) Indices

To calculate meteorological drought, in this study it was decided to choose the indices aSPI [16, 27] and RDI [16, 27, 44]. aSPI and RDI require only fP [14]) and rP–PET [18], respectively, as input data. In addition, these two indices of MeD are the most sensitive for semi-arid agricultural regions [16, 38, 48] such as the state of Sinaloa. In this study, the negative anomalies of aSPI and RDI were considered as dry events, as in Wang et al. [49].

Scales and Time Steps of Meteorological Drought (Med) Indices

To establish the scales and time steps for MeD, the associations between aSPI–RDI and the following were calculated: higher predictive sensitivity for rainfed maize yield [3 months (Nov–Jan) [16], semester with the highest annual contribution of rP [6 months (Jun–Nov) [26], quarter with the highest annual contribution of rP [3 months (Jul–Sep) [17]], semester with the highest annual contribution of PET [6 months (Mar–Aug) [17] and due to the follow-up order of agricultural years [16, 27, 50] [12 months (Oct–Sep).

Significant Trends of Meteorological Droughts (Med)

Significant Non-Parametric Trends (SNT)

To identify SNT (MeD > |1.96|), in each series of aSPI and RDI, the Mann–Kendall method was applied [51, 52]. To estimate the magnitudes of SNT, Sen’s slope method was applied [10; 51]. The annual results were multiplied by 10 to represent decadal changes.

Significant Linear Trends (SLT)

To reinforce the SNT study, a SLT analysis was also applied. This tool was based on the least squares method [53]. With the square root of the coefficient of determination (R2) of each SLT, Pe was determined [16]. Each was compared with Cpe = |0.329|, for n = 36 [54]. If Pe > |CPe|, then .

Climate Indices

Pacific Decadal Oscillation (PDO) and Atlantic Multidecadal Oscillation (AMO)

According to Ormaza et al. [55], PDO is the most important principal component of the monthly sea surface temperature anomalies in the North Pacific (20–65N and 120E–100W). The AMO is defined as the average of the sea surface temperature anomalies of the North Atlantic Ocean (0–70 °N) minus the average of the global sea surface temperature [56, 57]. PDO and AMO monthly data series were obtained from the National Oceanic and Atmospheric Administration [25], at https://psl.noaa.gov/ data/climateindices/list/.

Scales and Time Steps of the Pacific Decadal Oscillation (PDO) and Atlantic Multidecadal Oscillation (AMO)

To define the PDO and AMO indices, the time steps Mar–Aug (6-month) and Jul–Sep (3-month) were used. According to Llanes et al. [17], these two scales and time steps are the ones that register the highest annual accumulated contributions of PET and rP.

Statistical Pre-Treatment of Data

Prior to the principal component analysis and based on Llanes et al. [16] and Wüthrich and Merz [58], a Z standardized normalization was applied to all data series. Using the Shapiro–Wilk methods [59], Anderson–Darling [60], Lilliefors [61], and Jarque–Bera [62], the normality of all series was calculated.

Principal Component Analysis

To extract the most important indices, scales, and time steps, a principal component analysis was applied to all series with SNT (aSPI–3, Nov–Jan and aSPI–6, Jun– Nov) and SNT (aSPI–12, Oct–Sep; RDI–6, Mar–Aug and RDI–12, Oct–Sep): Z PC–1o was then extracted [6].

Association and Predictive Model Between the

First Principal Component (Z PC–1o) vs Pacific

Decadal Oscillation (Z PDO) and the Atlantic Multidecadal Oscillation (Z AMO)

To test the association of Z PC–1o vs Z PDO and Z AMO, a Pe correlation analysis was performed [63]. The correlations were Z PC–1o vs Z PDO–3 (Jul–Sep), Z PDO–6 (Mar–Aug), Z AMO–3 (Jul–Sep); and Z AMO–6 (Mar–Aug) was called Z PC–1c. The Pe values obtained were compared with Cpe = |0.329|, for n = 36 [53]. If Pe > |CPe|, Pe was considered significant. To predict Z PC–1o, a linear regression was applied, where Z PDO–3 (Jul–Sep), Z PDO–6 (Mar–Aug), Z AMO–3 (Jul–Sep) and Z AMO–6 (Mar –Aug) were the independent variables.

Validation of Model

A dispersion analysis was carried out between Z PC–1o vs Z PC–1c, and a Pe correlation was calculated. For the period 1981–2017, and after a Shapiro-Wilk normality analysis, linear regression of Z PC–1c was calculated. The coefficient of Pe of the series Z PC–1c was calculated by means of the square root of the coefficient of determination [Pe = (R2)0.5]. A hypothesis test was applied, where the Pe coefficients for Z PC–1o vs Z PC–1c, as well as for SLT of Z PC–1c, were compared with Cpe = |0.329|, for n = 36 [54]. If Pe > |CPe|, then Pe = significantly different from zero. This methodology was developed in accordance with what was described by Llanes et al. [16]. All statistical analyses were carried out with statistical significance = 0.05.

Software Used and Flowchart of the Applied Methodology

For this study, the programs used were DrinC 1.7 (91), XLstat 2023, PAST 4.08, Surfer 10 and CorelDRAW 2019.

In Fig. 2 (flowchart), each step carried out in the methodology of this study can be seen more clearly.

Fig. 2.

Flowchart of the methodology applied in this study.

https://www.pjoes.com/f/fulltexts/187162/PolJEnvironStud-2025-34-2-1515-1527_Fig_2_min.jpg

Results and Discussion

Significant Non-Parametric Trends (SNT)

Mann–Kendall Trend Identification

As shown in Table 1, three SNT were recorded, one negative [aSPI–3 (Nov–Jan) = −2.221 for Sanalona II] and two positive [aSPI–6 (Jun–Nov) = 2.003 for Potrerillos and RDI–6 (Mar–Aug) = 2.357 for Rosario]. The results from Table 1 can be attributed to the following annual average rP trends: a negative trend for Dec–Jan in Culiacán (20 km from Sanalona II) [64], a positive trend for Jun–Sep in Siqueros (17 km from Potrerillos) [64] and minimum accumulated rP with greater magnitude for Jan–Dec in Rosario (607.04 mm year−1) [17]. These results are in agreement with Sekhar et al. [65], Llanes et al. [66] and Nandi and Biswas [67], because for α = 0.05, significant non-parametric trends must comply with SNT > | 1.96|.

Table 1.

Identification of seasonal and annual significant non-parametric trends (SNT, framed magnitudes).

Weather stationaSPI–3 (Nov–Jan)aSPI–6 (Jun–Nov)aSPI–12 (Oct–Sep)RDI–3 (Jul–Sep)RDI–6 (Mar–Aug)RDI–12 (Oct–Sep)
Badiraguato–0.6000.6130.191–0.2040.899–0.341
Culiacán–0.9951.090–0.3001.3900.0950.436
El Playón–0.6690.490–0.4090.9130.7760.000
El Varejonal–0.926–0.9400.640–1.4851.022–0.313
Ixpalino–1.281–0.668–0.490–0.886–0.354–1.362
Jaina–1.648–0.8450.259–0.9541.3080.164
La Concha–1.7860.736–0.6130.804–0.681–0.300
Las Tortugas–0.9540.518–0.7630.586–0.708–0.518
Potrerillos0.8852.0031.0631.0900.2041.444
Rosario–0.7651.1170.4500.4092.3571.185
Sanalona II–2.221–0.749–1.498–0.490–0.913–1.731
Sta. Cruz de A.–0.709–1.2530.341–0.5041.4850.654

Magnitude of Significant Non–Parametric Trends (SNT) by Estimating Sen’s Slope

In Table 2, the SNT are presented for aSPI–3 (Nov– Jan) = −0.07 year−1 and aSPI–6 (Jun–Nov) = 0.05 year−1, for the Sanalona II and Potrerillos stations, respectively. These magnitudes from Table 2 are similar to those calculated by Syed et al. [6], who found an SNT for aSPI = −0.03 yr−1 for Saudi Arabia for the period 1985–2020, indicating a continuous increase in dry periods. According to Llanes et al. [16], these same dry periods that occurred in Sinaloa can considerably reduce crop yields, mainly rainfed [68] (such as corn) as well as aquifer recharge, as it is a region with a mostly semi-arid A.I. [69, 70].

Significant Linear Trends (SLT)

In Table 3 and for Sanalona II, the indices aSPI–3 (Nov–Jan), aSPI–12 (Oct–Sep) and RDI–12 (Oct–Sep), registered negative SLT (−0.0376, −0.0313 and −0.034, respectively), indicating that in this region (central Sinaloa) the planting–harvest dates of the crops for the autumn–winter and spring–summer cycles should be reassessed, since in future years yields may decrease significantly. RDI–3 (Jul–Sep) was the only index that did not register SLT (Table 3), which can be attributed to the fact that for Jul–Sep, the magnitude of the rainfall irregularity does not seriously affect rainfed crops [71]; that is, the trend of MeD is not significant [6; 72].

For Table 4, the most significant SLT was found for aSPI–3 (Nov–Jan) at the Sanalona II station (R2 = 0.169, r = 0.411), which indicates that for the central region of Sinaloa and for the period Nov –Jan, the greatest negative trends of MeD occurred [6; 16].

Normality values (p–values) for all data series were Shapiro Wilk (from 0.113 to 0.702), Anderson–Darling (from 0.158 to 0.824), Lilliefors (from 0.071 to 0.859) and Jarque–Bera (from 0.416 to 0.846).

Table 2.

Magnitude of seasonal and annual significant non-parametric trends (SNT, framed magnitudes).

Weather stationaSPI–3 (Nov–Jan)aSPI–6 (Jun–Nov)aSPI–12 (Oct–Sep)RDI–3 (Jul–Sep)RDI–6 (Mar–Aug)RDI–12 (Oct–Sep)
Badiraguato–0.0580.0780.113–0.0520.075–0.106
Culiacán–0.0730.060–0.0620.0590.0610.063
El Playón–0.0630.069–0.0340.0640.0700.043
El Varejonal–0.061–0.0500.044–0.0590.042–0.043
Ixpalino–0.057–0.072–0.088–0.052–0.092–0.082
Jaina–0.043–0.0540.051–0.0320.0540.042
La Concha–0.0700.083–0.1240.060–0.069–0.111
Las Tortugas–0.0440.046–0.0980.068–0.071–0.080
Potrerillos0.0460.0510.0720.0600.0950.089
Rosario–0.0550.0440.0580.0630.0540.056
Sanalona II–0.065–0.055–0.058–0.054–0.053–0.063
Sta. Cruz de A.–0.058–0.0470.079–0.0700.0910.090
Table 3.

Phase and magnitude of seasonal and annual significant linear trends (SLT, framed magnitudes).

Weather station (Variable y)Drought indices (Variable x)
aSPI–3 (Nov–Jan)aSPI–6 (Jun–Nov)aSPI–12 (Oct–Sep)RDI–3 (Jul–Sep)RDI–6 (Mar–Aug)RDI–12 (Oct–Sep)
Badiraguatoy = –0.010x y = 0.011x y = 0.011x y = –0.004x y = 0.020x y = 0.004x
Culiacány = –0.024x y = 0.012x y = –0.008x y = 0.020x y = 0.003x y = 0.002x
El Playóny = –0.015x y = 0.008x y = –0.004x y = 0.025x y = 0.012x y = 0.003x
El Varejonaly = –0.016x y = –0.020x y = 0.017x y = –0.026x y = 0.017x y = –0.002x
Ixpalinoy = –0.024x y = –0.012x y = –0.009x y = –0.017x y = 0.001x y = –0.014x
Jainay = –0.025x y = –0.013x y = 0.012x y = –0.014x y = 0.022x y = 0.006x
La Conchay = –0.027x y = 0.011x y = –0.013x y = 0.021x y = –0.008x y = –0.005x
Las Tortugasy = –0.014x y = 0.011x y = –0.018x y = 0.007x y = –0.019x y = –0.016x
Potrerillosy = 0.013x y = 0.036x y = 0.022x y = 0.027x y = 0.007x y = 0.024x
Rosarioy = –0.013x y = 0.011x y = 0.009x y = 0.017x y = 0.033x y = 0.022x
Sanalona IIy = –0.038x y = –0.014x y = –0.031x y = –0.017x y = –0.007x y = –0.034x
Sta. Cruz de A.y = –0.015x y = –0.018x y = 0.006x y = –0.012x y = 0.014x y = 0.006x
Table 4.

Coefficients of determination and correlation of seasonal and annual significant linear trends (SLT, framed magnitudes).

Drought indexBadiraguatoCuliacánEl PlayónEl VarejonalIxpalinoJainaLa ConchaLas TortugasPotrerillosRosarioSanalona IISta. Cruz de A.
aSPI–3 (Nov–Ene)R2 = 0.008R2 = 0.068R2 = 0.031R2 = 0.032R2 = 0.069R2 = 0.074R2 = 0.095R2 = 0.026R2 = 0.019R2 = 0.026R2 = 0.169R2 = 0.023
r = 0.089r = 0.261r = 0.176r = 0.179r = 0.263r = 0.272r = 0.308r = 0.161r = 0.138r = 0.161r = 0.411r = 0.152
aSPI–6 (Jun–Nov)R2 = 0.014R2 = 0.017R2 = 0.007R2 = 0.042R2 = 0.016R2 = 0.018R2 = 0.013R2 = 0.014R2 = 0.138R2 = 0.015R2 = 0.021R2 = 0.037
r = 0.118r = 0.130r = 0.084r = 0.205r = 0.126r = 0.134r = 0.114r = 0.118r = 0.371r = 0.122r = 0.145r = 0.192
aSPI–12 (Oct–Sep)R2 = 0.013R2 = 0.007R2 = 0.002R2 = 0.033R2 = 0.009R2 = 0.017R2 = 0.019R2 = 0.034R2 = 0.052R2 = 0.009R2 = 0.110R2 = 0.003
r = 0.114r = 0.084r = 0.045r = 0.182r = 0.095r = 0.130r = 0.138r = 0.184r = 0.228r = 0.095r = 0.332r = 0.055
RDI–3 (Jul–Sep)R2 = 0.002R2 = 0.043R2 = 0.068R2 = 0.074R2 = 0.032R2 = 0.021R2 = 0.045R2 = 0.005R2 = 0.079R2 = 0.030R2 = 0.031R2 = 0.015
r = 0.045r = 0.207r = 0.261r = 0.272r = 0.179r = 0.145r = 0.212r = 0.071r = 0.281r = 0.173r = 0.176r = 0.122
RDI–6 (Mar–Ago)R2 = 0.041R2 = 0.001R2 = 0.015R2 = 0.031R2 = 0.001R2 = 0.053R2 = 0.006R2 = 0.038R2 = 0.006R2 = 0.118R2 = 0.005R2 = 0.022
r = 0.202r = 0.032r = 0.122r = 0.176r = 0.032r = 0.230r = 0.077r = 0.195r = 0.077r = 0.344r = 0.071r = 0.148
RDI–12 (Oct–Sep)R2 = 0.002R2 = 0.001R2 = 0.001R2 = 0.001R2 = 0.020R2 = 0.004R2 = 0.002R2 = 0.027R2 = 0.063R2 = 0.054R2 = 0.125R2 = 0.004
r = 0.045r = 0.032r = 0.032r = 0.032r = 0.141r = 0.063r = 0.045r = 0.164r = 0.251r = 0.232r = 0.354r = 0.251
n = 36CPe =|0.329|

Principal Component Analysis

PC–1o (Equation 39) explained 47.99% of the variance of the drought trend. At the Sanalona II station, the indices with the greatest contribution to PC–1o (Equation 20) are presented. In descending order, they are Z aSPI–12 (Oct– Sep), Z RDI–12 (Oct–Sep) and Z aSPI–3 (Nov–Jan). The results from Equation 39 are in agreement with Llanes et al. [66], who calculated a SNT in the annual average rP for the Sanalona II station in the period 1982–2014.

(39)
PC–1o=0.26·[Z aSPI−3 (Jul– Sep, Sanalona II)]+0.02·[Z aSPI−6 (Jun– Nov, Potrerillos)]+0.37·[Z aSPI−12 (Oct– Sep, Sanalona II)]+0.02·[Z RDI–6 (Mar– Aug, Rosario)]+0.33·Z RDI−12 (Oct– Sep, Sanalona II)].

Association Between Z PC–1o vs Z Pacific Decadal

Oscillation (PDO) and Z Atlantic Multidecadal

Oscillation (AMO)

All variables [Z PC–1o vs Z PDO–3 (Jul–Sep), Z PDO–6 (Mar–Aug), Z AMO–3 (Jul–Sep) and Z AMO–6 (Mar–Aug)] presented normality, with Shapiro Wilk values (W) ranging from W = 0.951 to W = 0.979. In general, there were stronger associations between Z PC– 1o vs Z PDO than with respect to Z AMO. The highest correlation was obtained for Z PC–1o vs Z PDO–3 (Jul– Sep) = 0.390 (Table 5).

According to Table 5 the –Z PC–1o phase trend in Sinaloa is more influenced by –PDO phase anomalies (La Niña events) [23], than by +AMO phase anomalies (La Niña events) [17]. The time step Jul–Sep (period with the highest accumulated rP) is associated more with –Z PC– 1o phase anomalies than the time step Mar–Aug (period with the highest PET accumulated) [17].

These results are in agreement with Huang et al. [72] since these authors point out that +PDO anomalies (El Niño events) are highly related to the occurrence of intense hurricanes (absence of meteorological droughts) in the eastern Pacific Ocean (west of Sinaloa), which occur mainly in summer (Jul–Sep).

Model for the First Predicted Principal Component (Z PC–1c)

Following Table 5, the linear regression model for the calculation of Z PC–1c (Equation 40) is presented, where the explanatory variables are Z PDO–6 (Mar–Aug), Z AMO–6 (Mar–Aug) and Z AMO–3 (Jul–Sep):

(40)
Z PC−1c=0.01+0.34·[Z PDO−6(Mar−Aug)]−0.81·[Z AMO−6(Mar−Aug)]+0.50·[Z AMO−3(Jul−Sep)].

Testing the normality of the residuals of the model, we have W = 0.92, p = 0.815 and average = 0.0. The results from testing the normality of the residuals of the model (Equation 40) show that the data meet the normality assumptions for a linear regression [16].

Table 5.

Significant coefficients of Pearson correlation (Pe) (framed magnitudes) between Z PC–1o vs Z PDO and Z AMO.

Climate indexCoefficients of correlation (Pe)
Z PDO–6 (Mar–Aug)0.370
Z PDO–3 (Jul–Sep)0.390
Z AMO–6 (Mar–Aug)–0.378
Z AMO–3 (Jul–Sep)–0.318
n = 36; CPe = |0.329|; bold = significant coefficient

Validation of Model

The correlation coefficient between Z PC–1o vs Z PC–1c is Pe = 0.522, which was greater than CPe = ï0.329ï(Fig. 3a). The data series of Z PC–1c for the period 1981–2017 presented normality (W = 0.969, p–value = 0.394). The SLT of Z PC–1c was significant: R2 = 0.255 = Pe = 0.505 > CPe (Fig. 3b).

In Fig. 3a and in accordance with Llanes et al. [16] and OCR [54], and based on PDO–AMO indices (independent variables), the proposed model (–Z PC–1c, Equation 40) is sensitive (Pe > CPe, significantly different from zero) to predicting trends of MeD (aSPI–RDI indices, dependent variables) in the state of Sinaloa.

The periods 1997–1998, 2007–2008, 2009–2012 and 2016–2017 presented –Z PC–1c phase anomalies (Fig. 3b), which, according to Llanes et al. [16] and Llanes [17], are periods where negative anomalies (droughts) were also recorded in the following drought indices, scales and time steps: –aSPI–3 (Jul–Sep), –aSPI–6 (Mar–Aug), –RDI–3 (Jul–Sep), –RDI–6 (Mar–Aug) and –aSPI–3 (Nov–Jan) [15].

Conclusions

For the first time for Sinaloa and through linear regression, the MeD trend was sensitively modeled (aSPI–RDI, dependent variables), using PDO–AMO (Mar–Aug and Jul–Sep, independent variables). This study contributes by demonstrating the utility of a methodology that can be applied to any region of the world where limited meteorological data (temperature– precipitation) is available. One of the limitations of this study is that only the trends of the aSPI and RDI indices are addressed. Therefore, in future research, the use of other dependent variables is recommended; for example, standardized precipitation evapotranspiration index, rainfall deciles and percentiles, aridity anomaly index, crop moisture index, or hydrothermal coefficient, among other indexes of MeD. To try to increase the predictive sensitivity of the proposed model, it is also recommended to expand the study with other independent variables; for example, multivariate ENSO index, oceanic El Niño index, southern oscillation index, or North Atlantic oscillation, among other independent variables. In the validation, two Pe were significantly different from zero:

Fig. 3.

a) Dispersion between the first principal component observed (Z PC–1o) and calculated (Z PC–1c) and b) significant linear trend (SLT) of Z PC–1c.

https://www.pjoes.com/f/fulltexts/187162/PolJEnvironStud-2025-34-2-1515-1527_Fig_3_min.jpg

1) Pe of Z PC–1o vs Z PC–1c > CPe (Fig. 3a) and 2) Pe of SLT and for Z PC–1c > CPe (Fig. 3b). Based on the results obtained, an adjustment is recommended in the sowing dates in Sinaloa, which should avoid the months and the intense anomalies (> |1.0|) with phase +AMO (Jul–Sep) and –PDO (Mar–Aug), because they are generators of intense MeD. These results can help to schedule irrigation intelligently in one of the most important agricultural states in North America, called “the breadbasket of Mexico.”