Modeling the Rate of Change in Irrigation Water Quality Between Spring and Winter in Coastal Wells in the Municipality of Guasave, Sinaloa, Mexico

1. Introduction

Globally, groundwater is one of the most important freshwater resources [1]. Worldwide, approximately 20% of groundwater is used for agricultural irrigation [2,3]. As this water resource is essential, reliable, and constant [4], it supports the economic sector of any region [5]. In arid and semi-arid regions, the study of groundwater takes on greater importance due to its vulnerability [6,7]. This vulnerability is essentially due to the rate of change of physicochemical variables [8,9], which is accentuated when:

a)

Annual cumulative evapotranspiration values exceed precipitation values, especially in spring and winter [1,10].

b)

The cumulative precipitation of summer (July–September) is approximately three-quarters of the total annual accumulation [11].

c)

The irregularity of precipitation intensifies droughts [12], and consequently, groundwater becomes increasingly inefficient [13,14,15].

Groundwater inefficiency depends on multiple phenomena, such as recharge and dilution [4,7,14,15,16], salt dissolution and leaching [15,17], water–rock interaction [15,18], evaporation [19], and the rise and fall of the static water level (s.l.) [4,14]. According to [20,21], s.l. variation is highly related to extreme summer seasonal recharge [12] and is one of the factors most closely associated with variation in quality [9]. Irrigation water quality is evaluated, in accordance with [22,23], based on two parameters; electrical conductivity (E.C.) and sodicity. Multiple proposals have been developed to evaluate sodicity, three of which stand out: original sodium adsorption ratio (SAR_or), corrected (SAR_co), and adjusted sodium adsorption rate (SAR_ad) [5,7,24]. According to [24,25], these three formulations require the cations Na⁺, Ca²⁺ and Mg²⁺ for their calculation; however, only two of them account for anion contamination; CO₃^2– (SAR_co) and HCO₃^– (SAR_co and SAR_ad).

In Mexico, multiple efforts have also been made to evaluate groundwater quality for irrigation for the spring and winter seasons [24,26]. These two sampling seasons are essential for understanding, in the context of extreme summer precipitation, for example, the dissolution and leaching of salts derived from agrochemicals [2,11,13,15]. The semi-arid municipality of Guasave is located in the northern part of the state of Sinaloa [27,28]. Due to its high agricultural productivity, this municipality is also known as ‘‘the agricultural heart of Mexico’’ [29]. Although [19] evaluated irrigation water quality (E.C. vs. SAR_or) for coastal wells in Guasave, predictive models of the rate of quality change between spring and winter are still lacking.

In this study, two groundwater sampling campaigns were conducted for the coastal portion of the Guasave municipality: spring (June 2025) and winter (February 2026). In each campaign, 29 irrigation wells were sampled. For each sample, the following variables were measured: pH, E.C., total dissolved solids (TDS), temperature (Temp), salinity (Sal), s.l., Na⁺, Ca²⁺, Mg²⁺, K⁺, Cl^–, SO₄^2–, CO₃^2–, and HCO₃^–. SAR_or, SAR_co, and SAR_ad were calculated. For all variables, two ordinary and one repeated-measure [30] Pearson (Pr) and Spearman (Sr) correlation were applied. For spring and winter, a paired test [31] was applied and a Piper diagram [23] was constructed. To model quality for spring and winter, as well as the rate of quality change between the two seasons, multiple linear (MLR) and non-linear (MNR) regressions were applied, prioritizing efficiency and ease of application [32]. Observed and predicted quality were numerically classified using the nomogram from [22] and based on E.C. vs. each of SAR_or, SAR_co, and SAR_ad. To validate MLR and MNR, six and seven tests were applied, respectively.

The goal of this study is to model the rate of change between spring and winter of groundwater quality for irrigation in coastal wells in the municipality of Guasave, Sinaloa, Mexico.

The models in this study are easily to apply and can contribute to developing an efficient and effective method for preventing soil contamination [33] resulting from irrigation with saline and sodic groundwater. Preventing soil contamination can also foster optimal agricultural yields [18], helping to ensure food sovereignty and the socioeconomic development of one of the municipalities with the highest agricultural productivity in Mexico.

2. Materials and Methods

2.1. Study Area

2.1.1. Irrigated Agriculture, Geology, and Climatology

This study was conducted in the municipality of Guasave, which is considered “the agricultural heart of Mexico,” largely because approximately 181,542 hectares are irrigated; in other words, its economy relies to a large extent on activities associated with irrigated agriculture [29,34]. The crops primarily grown in this municipality are corn, potatoes, beans, green chili peppers, tomatoes, sorghum, and chickpeas [34].

The municipality of Guasave is composed of clastic sediments resting upon irregular igneous, sedimentary, and metamorphic rocks [19]. The deltaic plain is composed of sandy, silty, and clayey sediments, which originated during four epochs; the Oligocene, Miocene, Pliocene, and Pleistocene [35].

In the state of Sinaloa and for the year 2025, an average annual accumulated precipitation of 889.10 mm year⁻¹ was recorded, with the period June–October (the summer season) accounting for 97.64% of this precipitation [36]. The El Playón weather station (located approximately 47 km south of the city of Guasave) recorded an average annual accumulated precipitation of 463.03 mm year⁻¹ for the 1969–2018 period [28].

2.2. Location of the Study Area and Samples

This study focused on the coastal portion of the municipality of Guasave (Figure 1).

A water sample was randomly collected from each of 29 irrigation wells (Figure 1). Samples were collected from the same 29 wells during two seasons; spring (June 9–11, 2025) and winter (February 5–6, 2026), for a total of 58 samples.

2.3. Sample Storage

All samples were stored in 250 mL amber glass bottles at 4–10 °C [37]. Following measurement of s.l., to ensure the homogeneity of ionic concentration, each sample was collected after a minimum of three minutes of pumping. For metal measurement, each sample was acidified with 2.5 mL of 64.7% concentrated nitric acid [38].

2.4. Measurement and Calculation of Physicochemical Variables

Eight variables were measured in situ for each sample: TDS, Temp, pH, E.C., s.l., Na⁺, Sal, and Ca²⁺. In the laboratory, six additional variables were measured: Mg²⁺, K⁺, Cl⁻, SO₄^2–, CO₃^2– and HCO₃^–.

Three variables evaluating sodicity hazard were calculated; SAR_or (Equation 1), SAR_co (Equations 2 and 3), and SAR_ad (Equations 4 and 5). Equations 1 [24,39], 2 and 3 described by [40,41], and cited by [24], and 4 and 5 described by [25] were applied.

$${\mathrm{S}\mathrm{A}\mathrm{R}\mathrm{S}}_{\mathrm{o}\mathrm{r}}={\displaystyle \frac{{\mathrm{N}\mathrm{a}}^{+}}{\sqrt{\frac{{\mathrm{C}\mathrm{a}}^{2+}+{\mathrm{M}\mathrm{g}}^{2+}}{2}}}},$$

(1)

$${\mathrm{S}\mathrm{A}\mathrm{R}}_{\mathrm{c}\mathrm{o}}={\displaystyle \frac{{\mathrm{N}\mathrm{a}}^{+}}{\sqrt{\frac{{\mathrm{C}\mathrm{a}}^{2+}+{\mathrm{M}\mathrm{g}}^{2+}}{2}}}}⌊1.0+(8.4-{\mathrm{p}\mathrm{H}}_{\mathrm{c}})⌋,$$

(2)

$${\mathrm{p}\mathrm{H}}_{\mathrm{c}}=\left({\mathrm{p}\mathrm{k}\mathrm{\text{'}}}_2-{\mathrm{p}\mathrm{k}\mathrm{\text{'}}}_{\mathrm{c}}\right)+\mathrm{p}\mathrm{C}\mathrm{a}+\mathrm{p}\mathrm{A}\mathrm{l}\mathrm{k},$$

(3)

$${\mathrm{S}\mathrm{A}\mathrm{R}}_{\mathrm{a}\mathrm{d}}={\displaystyle \frac{{\mathrm{N}\mathrm{a}}^{+}}{\sqrt{\frac{{\mathrm{C}\mathrm{a}}_{\mathrm{e}\mathrm{q}}^{2+}+{\mathrm{M}\mathrm{g}}^{2+}}{2}}}},$$

(4)

$${\mathrm{C}\mathrm{a}}_{\mathrm{e}\mathrm{q}}^{2+}={10}^{\log \left(\mathrm{x}\right)}\times 0.17758,$$

(5)

2.5. Normality, Ordinary Correlation, Repeated Measures Correlation, and Correlation Hypothesis Testing for Physicochemical Variables

Prior to calculating correlations between physicochemical variables, a Shapiro–Wilk normality test was applied [31].

Pr was applied to the series that exhibited normality, while Sr was applied to those that did not. Pr and Sr were applied in three instances: a) ordinary correlation for spring, b) ordinary correlation for winter, and c) repeated measures correlation for spring and winter [30], which evaluated the seasonal influence of summer. The values involved in this last correlation were obtained through the following subtraction: values that were applied in repeated measures correlation = winter values – spring values.

Through hypothesis testing, significant correlations were identified that met the criteria of either Pr ≥ ?Pearson’s critical correlation (Pcr)? or Sr ≥ ?Spearman’s critical correlation (Scr)? [42]. The Pcr and Scr coefficients considered in this study were 1) for n = 23 (Pcr = 0.413 and Scr = 0.418) and 2) for n = 29 (Pcr = 0.367 and Scr = 0.370).

2.6. Paired Tests

Paired sample mean difference tests were applied [31]. If both series of values showed normality, the Student’s t-test was applied; if one or both series were non-normal, the Wilcoxon test was applied.

2.7. Hydrochemical Composition

To identify seasonal changes in the hydrochemical composition and characterization of the water, a Piper diagram was constructed for spring and for winter data. This diagram is one of the most widely used plots in the literature, as it is particularly useful for examining evolutions and ionic mixtures [23,43,44].

2.8. Observed Irrigation Water Quality Classification, and Prediction–Classification of Irrigation Water Quality

The Richards nomogram [22] was used to classify the observed (Figure 2a) and predicted qualities (Figure 2b). This classification evaluates quality using two parameters; E.C. and SAR [23,45].

Specifically for observed quality classification, an arbitrary numerical value was assigned to each, ranging from C1–S1 = 1 to C5–S5 = 25 (Figure 2a).

For the predicted quality, polynomial-type MLR models were used as the first option. When an MLR model failed to meet one of the validation criteria, a quadratic polynomial-type MNR model was used as the second option. For the MLR and MNR dependent variables, the observed numerical values for three water quality parameters were used; E.C. vs. SAR_or, E.C. vs. SAR_co, and E.C. vs. SAR_ad. Prioritizing ease of application, it was decided to use four and six independent variables; pH, E.C., Temp, and s.l. for MLR and pH, E.C., TDS, Temp, Sal, and s.l. for MNR. Once the numerical values for predicted quality were obtained, they were also classified using equivalence intervals from C1–S1 = 0.01–1.00 to C5–S5 = 24.01–25.00 (Figure 2b). Quality was predicted for three conditions: a) spring, b) winter, and c) difference between winter and spring.

2.9. Salinity (Sal) Estimation Regressions

It was decided to provide predictive models for Sal, based on E.C. This is in the case where salt meters in the wells are dispensed with.

To predict Sal, simple linear regressions were first applied. When normality was not observed in the residuals of these linear regressions, simple non-linear regressions were used as a second option, through a quadratic polynomial. Since the observed and predicted salinity series did not exhibit normality, Sr was applied. For the predicted Sal models to be considered significant, the condition Sr ≥ ?Scr? had to be met (Scr = 0.418, n = 23).

2.10. Validation of Predictive Quality Models

2.10.1. Normality, Homoscedasticity, Homogeneity and Independence of Residuals, Minimal Multicollinearity and Performance Metrics

Four validation tests were applied to the MLR residuals: 1) Shapiro–Wilk normality [46], 2) homoscedasticity (Breusch–Pagan and White) [47,48], 3) homogeneity [49], and 4) independence (Durbin–Watson) [50], (Box–Pierce) [51], (Ljung–Box) [51], and (McLeod–Li) [52].

Furthermore, for the MLR, it was verified that all independent variables exhibited minimal multicollinearity—specifically, a low variance inflation factor (VIF < 3.0) [53,54].

Seven performance tests were applied to the MLR and MNR; Sr (fit), coefficient of determination (${\mathrm{R}}^2={\mathrm{S}\mathrm{r}}^2$), root mean square error (RMSE), mean square error (MSE), Nash–Sutcliffe efficiency (NSE), percent bias (PBIAS), and mean error (ME). These tests are widely used to evaluate the performance of predictive models [55,56,57,58,59,60]. According to [42,61], the optimal performances are Sr $\ge$ ?Scr? (Scr = 0.418, n = 23), R² $\ge$ ${0.418}^2$, RMSE $~$ 0.000, MSE $~$ 0.000, 0.750 $<$ NSE $\le$ 1.000, PBIAS $<$ $\pm$ 10.000, and ME $~$ 0.000.

2.11. Software, Equipment, and Statistical Significance

The following software programs were used in this study: Paleontological Statistics (PAST) version 5.3, Microsoft Office 365, XLSTAT Advanced version 2026.1.0, and CorelDRAW 2019 version 21.0.0.593.

The following equipment was also used: a Garmin Map 62s GPS, two Hanna Instruments photometers models HI83399–01 (for Mg²⁺, K⁺, and Cl⁻), and HI83325–01 (for SO₄²⁻), a Hanna Instruments potentiometer HI98130 (for E.C., pH, Temp, and TDS), three Horiba LAQUAtwin ion-selective meters (for Na₊, Ca²⁺, and Sal), and a level probe and Truper tape measure. CO₃²⁻ and HCO₃⁻ were measured via volumetry, in accordance with [62].

The significance level for all statistical analyses isα = 0.05.

3. Results

3.1. Spring

3.1.1. Spatial Variability of Physicochemical Parameters

As shown in Table 1, the maximum, minimum, and average values for the physicochemical parameters are 8.23, 6.90, and 7.59 (pH); 14,910.00 μS cm⁻¹, 560.00 μS cm⁻¹, and 2,193.45 μS cm⁻¹ (E.C.); 7.46 ppt, 0.27 ppt, and 1.09 ppt (TDS); 40.10 °C, 27.00 °C, and 30.03 °C (Temp); 8.20% NaCl, 0.30% NaCl, and 1.26% NaCl (Sal); 5.89 m, 1.98 m, and 3.42 m (s.l.); and 117.44 mEq L⁻¹, 2.26 mEq L⁻¹, and 23.54 mEq L⁻¹ (Na⁺). The further maximum, minimum and average values are 20.96 mEq L^–1, 0.42 mEq L^–1, and 6.31 mEq L^–1 (Ca²⁺); 33.74 mEq L^–1, 0.66 mEq L^–1, and 3.03 mEq L^–1 (Mg²⁺); 1.38 mEq L^–1, 0.04 mEq L^–1, and 0.15 mEq L^–1 (K⁺); 22.99 mEq L^–1, 0.42 mEq L^–1 and 3.42 mEq L^–1 (Cl^–); 14.16 mEq L^–1, 0.35 mEq L^–1, and 4.39 mEq L^–1 (SO₄^2–); all values zero for CO₃^2–; 13.33 mEq L^–1, 6.66 mEq L^–1, and 10.25 mEq L^–1 (HCO₃^–); 31.25, 1.33, and 12.55 (SAR_or ); 74.11, 3.29, and 31.00 (SAR_co); and 37.53, 2.09, and 17.03 (SAR_ad) (Table 1).

3.1.2. Observed Irrigation Water Quality Classification

Based on the E.C. vs. SAR_or results (Figure 3a), the following quality classifications were recorded: one well with C5–S3 (well 20), one well with C3–S5 (well 17), six wells with C4–S4 (wells 1, 6, 9–11, and 21), four wells with C3–S4 (wells 13, 18, 28, and 29), two wells with C3–S3 (wells 12 and 24), eight wells with C3–S2 (wells 5, 8, 14–16, 22, 23, and 25), four wells with C3–S1 (wells 4, 7, 19, and 27), and three wells with C2–S1 (wells 2, 3, and 26). Figure 3b (E.C. vs. SAR_co) shows the following quality classifications: one well with C5–S5 (well 20), six wells with C4–S5 (wells 1, 6, 9–11, and 21), six wells with C3–S5 (wells 13, 17, 18, 24, 28, and 29), eight wells with C3–S4 (wells 5, 8, 12, 14, 16, 22, 23, and 25), two wells with C3–S3 (wells 15 and 27), three wells with C3–S2 (wells 4, 7, and 19), two wells with C2–S2 (wells 2 and 3), and one well with C2–S1 (well 26). Figure 3c (E.C. vs. SAR_ad) shows the following quality classifications: one well with C5–S4 (well 20), one well with C4–S5 (well 9), two wells with C3–S5 (wells 17 and 28), five wells with C4–S4 (wells 1, 6, 10, 11, and 21), five wells with C3–S4 (wells 12, 13, 18, 24, and 29), four wells with C3–S3 (wells 5, 22, 23, and 25), seven wells with C3–S2 (wells 7, 8, 14–16, 19, and 27), one well with C3–S1 (well 4), and three wells with C2–S1 (wells 2, 3, and 26).

3.2. Winter

3.2.1. Spatial Variability of Physicochemical Parameters

As shown in Table 2, the maximum, minimum, and average values for the parameters are 8.09, 6.64, and 7.18 (pH); 10,230.00 μS cm⁻¹, 600.00 μS cm⁻¹; and 1,849.31 μS cm⁻¹ (E.C.); 5.13 ppt, 0.31 ppt, and 0.93 ppt (TDS); 29.10 °C, 21.70 °C, and 26.71 °C (Temp); 8.30% NaCl, 0.40% NaCl, and 1.22% NaCl (Sal); 6.80 m, 1.29 m, and 2.30 m (s.l.); and 213.14 mEq L⁻¹, 9.57 mEq L⁻¹, and 38.08 mEq L⁻¹ (Na⁺). The maximum, minimum, and average values for the remaining parameters are 23.95 mEq L^–1, 0.55 mEq L⁻¹, and 6.72 mEq L⁻¹ (Ca²⁺); 40.32 mEq L^–1, 0.74 mEq L^–1, and 5.00 mEq L^–1 (Mg²⁺); 1.85 mEq L^–1, 0.04 mEq L⁻¹, and 0.19 mEq L^–1 (K⁺); 119.88 mEq L^–1, 0.04 mEq L⁻¹, and 7.13 mEq L^–1 (Cl^–); 13.12 mEq L^–1, 0.92 mEq L^–1, and 3.66 mEq L^–1 (SO₄^2–); all values zero for CO₃^2–; 51.74 mEq L⁻¹, 7.45 mEq L^–1, and 17.83 mEq L^–1 (HCO₃^–); 52.07, 3.60, and 18.51 (SAR_or); 89.92, 3.99, and 34.06 (SAR_co); and 66.91, 6.24, and 24.62 (SAR_ad) (Table 2).

3.2.2. Observed Irrigation Water Quality Classification

The quality classifications for E.C. vs. SAR_or (Figure 4a) are: one well with C5–S5 (well 20), four wells with C3–S5 (wells 9, 13, 18, and 28), four wells with C4–S4 (wells 1, 6, 21, and 29), four wells with C3–S4 (wells 5, 11, 12, and 24), seven wells with C3–S3 (wells 8, 10, 14, 22, 23, 25, and 27), seven wells with C3–S2 (wells 2–4, 7, 16, 17, and 19), one well with C3–S1 (well 26), and one well with C2–S1 (well 15). For E.C. vs. SAR_co (Figure 4b), the quality classifications are: one well with C5–S5 (well 20), four wells with C4–S5 (wells 1, 6, 21, and 29), nine wells with C3–S5 (wells 9–13, 17, 18, 24, and 28), seven wells with C3–S4 (wells 5, 8, 14, 16, 22, 23, and 25), one well with C3–S3 (well 27), five wells with C3–S2 (wells 2–4, 7, and 19), one well with C3–S1 (well 26), and one well with C2–S3 (well 15). Figure 4c (E.C. vs. SAR_ad) shows the following quality classifications: one well with C5–S5 (well 20), three wells with C4–S5 (wells 6, 21, and 29), four wells with C3–S5 (wells 9, 13, 18, and 28), one well with C4–S4 (well 1), seven wells with C3–S4 (wells 5, 10–12, and 23–25), eight wells with C3–S3 (wells 2, 3, 8, 14, 16, 17, 22, and 27), four wells with C3–S2 (wells 4, 7, 19, and 26), and one well with C2–S2 (well 15).

3.3. Correlational Analysis Between Physicochemical Variables

3.3.1. Ordinary Correlation for Spring

As presented in Table 3, the following Sr values for E.C. were significant: 0.997 (TDS), 0.933 (Sal), –0.463 (s.l.), 0.859 (Na⁺), 0.426 (Ca²⁺), 0.544 (Mg²⁺), 0.775 (Cl^–), 0.843 (SO₄²⁻), 0.619 (SAR_or), 0.695 (SAR_co), and 0.668 (SAR_ad).

SAR_or showed the following significant correlations: Pr = 0.489 (pH), Sr = 0.619 (E.C.), Sr = 0.611 (TDS), Sr = 0.531 (Sal), Sr = 0.880 (Na⁺), Sr = 0.563 (SO₄^2–), and Pr = 0.541 (HCO₃⁻). SAR_co was significantly correlated with E.C. (Sr = 0.695), TDS (Sr = 0.687), Sal (Sr = 0.604), Na⁺ (Sr = 0.927), Cl^– (Sr = 0.406), SO₄^2– (Sr = 0.638), and HCO₃^– (Pr = 0.489). SAR_ad was significantly correlated with pH (Pr = 0.452), E.C. (Sr = 0.668), TDS (Sr = 0.659), Sal (Sr = 0.576), Na⁺ (Sr = 0.909), Cl^– (Sr = 0.398), SO₄^2– (Sr = 0.610), and HCO₃^– (Pr = 0.545).

Temp and K⁺ showed no significant correlation with E.C., SAR_or, RAS_co, or SAR_ad.

3.3.2. Ordinary Correlation for Winter

As shown in Table 4, the significant values of Sr and for E.C. are 0.999 (TDS), 0.992 (Sal), 0.738 (Na⁺), 0.482 (Mg²⁺), 0.906 (Cl^–),0.635 (SO₄^2–),0.377 (SAR_or), and 0.499 (SAR_co).

SAR_or showed the following significant correlations for Sr: 0.741 (pH), 0.377 (E.C.), 0.377 (TDS), 0.860 (Na⁺), –0.625 (Ca²⁺), and –0.389 (Mg²⁺). SAR_co showed the following significant correlations: Pr = 0.708 (pH), Sr = 0.499 (E.C.), Sr = 0.501 (TDS), Sr = 0.479 (Sal), Sr = 0.867 (Na⁺), Sr = −0.370 (Ca²⁺), Sr = 0.424 (Cl^–), and Sr = 0.511 (SO₄^2–). The significant correlations for Sr and for SAR_ad are 0.707 (pH), 0.370 (TDS), 0.826 (Na⁺), −0.576 (Ca²⁺), and –0.447 (Mg²⁺).

Temp, s.l., K⁺, and HCO₃^– showed no significant correlations with E.C., SAR_or, SAR_co, or SAR_ad.

3.3.3. Repeated Measures Correlation (Difference Between Winter and Spring)

Table 5 shows the following significant correlations with E.C.: TDS (0.977), Sal (0.563), and SO₄²⁻ (0.773).

As seen in Table 5, the significant Sr correlations are 0.671 for SAR_or vs. Na⁺ and −0.393 for SAR_or vs. Mg²⁺, 0.524 for SAR_co vs. Mg²⁺ and 0.605 for SAR_co vs. HCO₃^–, and 0.623 for SAR_ad vs. Na⁺ and −0.455 for SAR_ad vs. Mg²⁺.

S.l., Temp, K⁺, and SO₄²⁻ showed no significant correlation with E.C., SAR_or, SAR_co, or SAR_ad.

3.4. Paired Tests

Of 17 physicochemical variables (Table 6), only four did not show significant changes in the means: Sal (p-value = 0.101), Ca²⁺ (p-value = 0.503), Cl^– (p-value = 0.756), and CO₃^2– (p-value = identical samples).

3.5. Hydrochemical Composition

3.5.1. Spring

As seen in Figure 5, out of 29 samples, 18 have a sodium-bicarbonate composition, 8 are sodium chloride-sulfate, 2 are magnesium-calcium-sulfate-chloride, and 1 is a magnesium-calcium-bicarbonate. All samples were associated with two cations: Na⁺ (26 samples) and Ca²⁺ (3 samples). Twenty-three associations with anions were observed: 19 with HCO₃^– and 4 with Cl^–, while 6 samples showed no dominant anion.

3.5.2. Winter

As shown in the Piper diagram in Figure 6, 22 wells had sodium-bicarbonate water, 5 had sodium-chloride-sulfate water, and 2 had calcium-magnesium-bicarbonate water. Twenty-seven samples were associated with Na⁺, and two samples showed no dominant cations. There were 24 associations with HCO₃^– and two with Cl^–. Three samples showed no dominant anion.

3.6. Quality Predictive Models

Nine models were obtained for quality prediction: four MLR models (Equations 6 and 8–10) and five MNR models (Equations 7 and 11–14).

$$\begin{gathered} \mathrm{S}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{q}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}(\mathrm{E}.\mathrm{C}.\mathrm{v}\mathrm{s}{\mathrm{S}\mathrm{A}\mathrm{R}}_{\mathrm{o}\mathrm{r}})=10.636+\left(1.084·\mathrm{p}\mathrm{H}\right)+\left(8.528\times {10}^{-4}·\mathrm{E}.\mathrm{C}.\right) \\ -\left(8.473\times {10}^{-2}·\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}\right)-\left(1.100·\mathrm{s}.\mathrm{l}.\right), \end{gathered}$$

(6)

$$\begin{gathered} \mathrm{S}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{q}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}\left(\mathrm{E}.\mathrm{C}.\mathrm{v}\mathrm{s}{\mathrm{S}\mathrm{A}\mathrm{R}}_{\mathrm{c}\mathrm{o}}\right)=-41.560+\left(16.798·\mathrm{p}\mathrm{H}\right)+\left(0.021·\mathrm{E}.\mathrm{C}.\right)-\left(34.443·\mathrm{T}\mathrm{D}\mathrm{S}\right)-\left(0.645·\mathrm{T}\mathrm{e}\mathrm{m}\right)+ \\ \left(0.057·\mathrm{S}\mathrm{a}\mathrm{l}\right)-\left(2.973·\mathrm{s}.\mathrm{l}.\right)-\left(1.010·{\mathrm{p}\mathrm{H}}^2\right)-\left(2.344\times {10}^{-6}·{\mathrm{E}.\mathrm{C}.}^2\right)+\left(9.458·{\mathrm{T}\mathrm{D}\mathrm{S}}^2\right)+\left(0.009·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}^2\right)- \\ \left(0.607·{\mathrm{S}\mathrm{a}\mathrm{l}}^2\right)+\left(0.315·{\mathrm{s}.\mathrm{l}.}^2\right), \end{gathered}$$

(7)

$$\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{q}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}(\mathrm{E}.\mathrm{C}.\mathrm{v}\mathrm{s}{\mathrm{S}\mathrm{A}\mathrm{R}}_{\mathrm{a}\mathrm{d}})=9.264+\left(1.361·\mathrm{p}\mathrm{H}\right)+\left(8.285\times {10}^{-4}·\mathrm{E}.\mathrm{C}.\right)-\left(0.102·\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}\right)-\left(1.046·\mathrm{s}.\mathrm{l}.\right),$$

(8)

$$\mathrm{W}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r} \mathrm{q}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}2026\left(\mathrm{E}.\mathrm{C}.\mathrm{v}\mathrm{s}{\mathrm{S}\mathrm{A}\mathrm{R}}_{\mathrm{o}\mathrm{r}}\right)=-16.273+\left(2.728·\mathrm{p}\mathrm{H}\right)+\left(1.519\times {10}^{-3}·\mathrm{E}.\mathrm{C}.\right)+\left(0.307·\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}\right)+\left(1.191·\mathrm{s}.\mathrm{l}.\right),$$

(9)

$$\mathrm{W}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r} \mathrm{q}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}\left(\mathrm{E}.\mathrm{C}.\mathrm{v}\mathrm{s}{\mathrm{S}\mathrm{A}\mathrm{R}}_{\mathrm{c}\mathrm{o}}\right)=-6.312+\left(2.250·\mathrm{p}\mathrm{H}\right)+\left(1.453\times {10}^{-3}·\mathrm{E}.\mathrm{C}.\right)+\left(0.104·\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}\right)-\left(0.076·\mathrm{s}.\mathrm{l}.\right),$$

(10)

$$\begin{gathered} \mathrm{W}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r} \mathrm{q}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}\left(\mathrm{E}.\mathrm{C}.\mathrm{v}\mathrm{s}{\mathrm{S}\mathrm{A}\mathrm{R}}_{\mathrm{a}\mathrm{d}}\right)=-183.166+\left(31.649·\mathrm{p}\mathrm{H}\right)+\left(0.077·\mathrm{E}.\mathrm{C}.\right)-\left(146.068·\mathrm{T}\mathrm{D}\mathrm{S}\right)+\left(5.111·\mathrm{T}\mathrm{e}\mathrm{m}\right)- \\ \left(1.799·\mathrm{S}\mathrm{a}\mathrm{l}\right)+\left(1.642·\mathrm{s}.\mathrm{l}.\right)-\left(1.976·{\mathrm{p}\mathrm{H}}^2\right)-\left(2.654\times {10}^{-5}·{\mathrm{E}.\mathrm{C}.}^2\right)+\left(108.782·{\mathrm{T}\mathrm{D}\mathrm{S}}^2\right)-\left(0.101·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}^2\right)- \\ \left(1.360·{\mathrm{S}\mathrm{a}\mathrm{l}}^2\right)-\left(0.175·{\mathrm{s}.\mathrm{l}.}^2\right), \end{gathered}$$

(11)

$$\begin{gathered} \mathrm{R}\mathrm{a}\mathrm{t}\mathrm{e} \mathrm{o}\mathrm{f} \mathrm{c}\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\left(\mathrm{E}.\mathrm{C}.\mathrm{v}\mathrm{s}{\mathrm{S}\mathrm{A}\mathrm{R}}_{\mathrm{o}\mathrm{r}}\right)=-242.361+(159.494·{\mathrm{p}\mathrm{H}}_{\mathrm{s}\mathrm{p}})-(0.009·{\mathrm{E}.\mathrm{C}.}_{\mathrm{s}\mathrm{p}})+(7.533·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{s}\mathrm{p}})- \\ (16.115·{\mathrm{s}.\mathrm{l}.}_{\mathrm{s}\mathrm{p}})-(183.746·{\mathrm{p}\mathrm{H}}_{\mathrm{w}\mathrm{i}})+(0.013·{\mathrm{E}.\mathrm{C}.}_{\mathrm{w}\mathrm{i}})+(17.126·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{w}\mathrm{i}})+(7.4747·{\mathrm{s}.\mathrm{l}.}_{\mathrm{w}\mathrm{i}})- \\ (10.277·{\mathrm{p}\mathrm{H}}_{\mathrm{s}\mathrm{p}}^2)+(9.048\times {10}^{-7}·{\mathrm{E}.\mathrm{C}.}_{\mathrm{s}\mathrm{p}}^2)-(0.126·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{s}\mathrm{p}}^2)+(2.250·{\mathrm{s}.\mathrm{l}.}_{\mathrm{s}\mathrm{p}}^2)+(12.051·{\mathrm{p}\mathrm{H}}_{\mathrm{w}\mathrm{i}}^2)- \\ (1.957\times {10}^{-6}·{\mathrm{E}.\mathrm{C}.}_{\mathrm{w}\mathrm{i}}^2)-(0.325·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{w}\mathrm{i}}^2)-(1.095·{\mathrm{s}.\mathrm{l}.}_{\mathrm{w}\mathrm{i}}^2), \end{gathered}$$

(12)

$$\begin{gathered} \mathrm{R}\mathrm{a}\mathrm{t}\mathrm{e} \mathrm{o}\mathrm{f} \mathrm{c}\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\left(\mathrm{E}.\mathrm{C}.\mathrm{v}\mathrm{s}{\mathrm{S}\mathrm{A}\mathrm{R}}_{\mathrm{c}\mathrm{o}}\right)=-30.007+(34.212·{\mathrm{p}\mathrm{H}}_{\mathrm{s}\mathrm{p}})-(0.008·{\mathrm{E}.\mathrm{C}.}_{\mathrm{s}\mathrm{p}})+(7.501·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{s}\mathrm{p}})- \\ (19.094·{\mathrm{s}.\mathrm{l}.}_{\mathrm{s}\mathrm{p}})-(77.177·{\mathrm{p}\mathrm{H}}_{\mathrm{w}\mathrm{i}})+(0.013·{\mathrm{E}.\mathrm{C}.}_{\mathrm{w}\mathrm{i}})+(8.053·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{w}\mathrm{i}})+(7.979·{\mathrm{s}.\mathrm{l}.}_{\mathrm{w}\mathrm{i}})-(2.380·{\mathrm{p}\mathrm{H}}_{\mathrm{s}\mathrm{p}}^2)+ \\ (8.264\times {10}^{-7}·{\mathrm{E}.\mathrm{C}.}_{\mathrm{s}\mathrm{p}}^2)-(0.120·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{s}\mathrm{p}}^2)+(2.600·{\mathrm{s}.\mathrm{l}.}_{\mathrm{s}\mathrm{p}}^2)+(5.124·{\mathrm{p}\mathrm{H}}_{\mathrm{w}\mathrm{i}}^2)-(1.902\times {10}^{-6}·{\mathrm{E}.\mathrm{C}.}_{\mathrm{w}\mathrm{i}}^2)- \\ (0.162·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{w}\mathrm{i}}^2)-(1.173·{\mathrm{s}.\mathrm{l}.}_{\mathrm{w}\mathrm{i}}^2), \end{gathered}$$

(13)

$$\begin{gathered} \mathrm{R}\mathrm{a}\mathrm{t}\mathrm{e} \mathrm{o}\mathrm{f} \mathrm{c}\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\left(\mathrm{E}.\mathrm{C}.\mathrm{v}\mathrm{s}{\mathrm{S}\mathrm{A}\mathrm{R}}_{\mathrm{a}\mathrm{d}}\right)=-218.100+(158.615·{\mathrm{p}\mathrm{H}}_{\mathrm{s}\mathrm{p}})-(0.008·{\mathrm{E}.\mathrm{C}.}_{\mathrm{s}\mathrm{p}})+(6.192·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{s}\mathrm{p}})- \\ (16.131·{\mathrm{s}.\mathrm{l}.}_{\mathrm{s}\mathrm{p}})-(159.999·{\mathrm{p}\mathrm{H}}_{\mathrm{w}\mathrm{i}})+(0.013·{\mathrm{E}.\mathrm{C}.}_{\mathrm{w}\mathrm{i}})+(10.582·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{w}\mathrm{i}})+(9.430·{\mathrm{s}.\mathrm{l}.}_{\mathrm{w}\mathrm{i}})- \\ (10.396·{\mathrm{p}\mathrm{H}}_{\mathrm{s}\mathrm{p}}^2)+(8.524\times {10}^{-7}{·\mathrm{E}.\mathrm{C}.}_{\mathrm{s}\mathrm{p}}^2)-(0.105·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{s}\mathrm{p}}^2)+(2.144·{\mathrm{s}.\mathrm{l}.}_{\mathrm{s}\mathrm{p}}^2)+(10.597·{\mathrm{p}\mathrm{H}}_{\mathrm{w}\mathrm{i}}^2)- \\ (1.903\times {10}^{-6}{·\mathrm{E}.\mathrm{C}.}_{\mathrm{w}\mathrm{i}}^2)-(0.204·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{w}\mathrm{i}}^2)-(1.257·{\mathrm{s}.\mathrm{l}.}_{\mathrm{w}\mathrm{i}}^2), \end{gathered}$$

(14)

3.6.1. Salinity (Sal) Estimation Regressions

To predict Sal content based on E.C., two simple non-linear regressions were obtained (Equations 15 and 16). For observed and predicted spring Sal (Equation 15), the p-values were p(normal) = 7.903$\times {10}^{-8}$, p(normal) = 6.496$\times {10}^{-8}$, and Sr = 0.905, respectively. For observed and predicted salinity in winter (Equation 16), the p-values were p(normal) = 4.583$\times {10}^{-8}$, p(normal) = 3.658$\times {10}^{-8}$, and Sr = 0.993, respectively.

$$\mathrm{S}\mathrm{a}\mathrm{l}\left(2025\right)=0.086+\left(5.367\times {10}^{-4}·\mathrm{E}.\mathrm{C}.\right)+(5.676\times {10}^{-10}·{\mathrm{E}.\mathrm{C}.}^2),$$

(15)

$$\mathrm{S}\mathrm{a}\mathrm{l}\left(2026\right)=0.079+\left(5.206\times {10}^{-4}·\mathrm{E}.\mathrm{C}.\right)+(2.753\times {10}^{-8}·{\mathrm{E}.\mathrm{C}.}^2),$$

(16)

3.7. Validation of Predictive Quality Models

3.7.1. Normality, Homoscedasticity, Homogeneity and Independence of Residuals, Minimal Multicollinearity and Performance Metrics

In the results shown in Table 7, the residuals exhibiting normality ranged from p(normal) = 0.063 [winter (E.C. vs. SAR_or)] to p(normal) = 0.274 [winter (E.C. vs. SAR_co)]. Two models were non-normal; p(normal) = 0.030 [winter (E.C. vs. SAR_ad)] and p(normal) = 0.002 [winter–spring (E.C. vs. SAR_or)].

Homoscedasticity ranged from p-value = 0.079 [Breusch–Pagan test for spring (E.C. vs. SAR_or)] to p-value = 0.878 [Breusch–Pagan test for winter (E.C. vs. SAR_co)]. Only one series was heteroscedastic: p-value = 0.040 (Breusch–Pagan test for spring [E.C. vs. SAR_co]).

All series that exhibited normality and homoscedasticity also exhibited homogeneity; that is, the means of the residuals were 0.000.

Models with independence ranged from p(value) = 0.117 [Ljung–Box test for winter–spring (E.C. vs SAR_co)] to p(value) = 0.942 [Durbin–Watson test for spring (E.C. vs SAR_or)]. One series showed autocorrelation, with a p-value of 0.036 [McLeod–Li test and for winter–spring (E.C. vs. SAR_ad)] (Table 7).

3.7.2. Minimal Multicollinearity

As shown in in Table 7, the minimal multicollinearity of the predictors ranged from VIF = 1.053 (for Temp) to VIF = 1.372 (for E.C.) across two models; [winter (E.C. vs. SAR_or) and spring (E.C. vs. SAR_ad)]. One series showed severe multicollinearity; namely VIF = 25.412 (for E.C.) in [winter–spring (E.C. vs. SAR_co)].

3.7.3. Performance Metrics

As seen in Table 7 and Figure 7, the correlations ranged from Sr = 0.470 (R² = 0.221) for winter–spring (E.C. vs. SAR_ad, Figure 7i) to Sr = 0.922 (R² = 0.851) for spring (E.C. vs. SAR_co, Figure 7b). Further performance metrics, the maximum and minimum values, respectively, were 2.033 (Figure 7a) and 0.804 (Figure 7b) for RMSE, 4.134 (Figure 7a) and 0.646 (Figure 7b) for MEC, 1.000 (Figure 7g and Figure 7i.991 (Figure 7d) for NSE, 4.628 $\times {10}^{-6}$ (Figure 7b) and 2.672$\times {10}^{-14}$ (Figure 7a) for PBIAS, and 7.405$\times {10}^{-7}$ (Figure 7b) and 3.669$\times {10}^{-16}$ (Figure 7i) for ME.

4. Discussion

According to [63], normal and non-alkaline water values (6.50 < pH < 8.40) like those shown in Table 1 and Table 2 indicate good nutrient availability [24]. These values indicating the lack of alkalinity are also decisive for the occurrence of zero CO₃^2– concentrations [63]. In general, the decreases in the values of E.C., TDS, Temp, Sal, and SO₄²⁻ (Table 2) can be attributed to the fact that summer precipitation, accounting for 72.96% of the total annual accumulation [11], tends to recharge the Sinaloa river aquifer [19,35]. This recharge dilutes the ionic concentration of the water [19] and decreases the depth of s.l. [35]. Table 2 presents, in general terms, an incremental effect on the average values of six ions (Na⁺, Ca²⁺, Mg²⁺, K⁺, Cl^–, and HCO₃^–) and three indices (SAR_or, SAR_co, and SAR_ad). This incremental effect can be attributed to heterogeneous hydrochemical responses such as mineral dissolution [24,64], soil salt leaching [65], cation exchange [65], and mixing with seawater [66].

In this study, 43 samples (73.13% of the total; Figure 3 and Figure 4) align with the results found by [23], as these authors also identified a high salinity hazard (C3 classification). Waters classified as C3 can only be used to irrigate alfalfa, wheat, oats, and barley [24], as these are salt-tolerant rather than salt-sensitive crops [67]. Based on Figure 4, to avoid water with a high sodification hazard (S4 classification) and thereby prevent a negative impact on soil structure [68] and undesirable crop development [69], it is suggested to use the water during winter, as this tends to lower SAR_or, SAR_co, and SAR_ad values.

As shown in Table 3, the highest correlations for E.C. were obtained with s.l., Na⁺, Ca²⁺, Mg²⁺, and SO₄^2–, and for SAR with Na⁺, Cl^–, SO₄^2–, and HCO₃^–. These correlations can be attributed to an increase in ionic concentration due to evaporative and minimal recharge effects [17,19].

According to [17], the highest correlations in Table 4 [(E.C. vs. Cl^–) and (SAR vs. pH, Ca²⁺, and Mg²⁺)] can be attributed to the effects of soil leaching, cation exchange, and salt dissolution. According to [70], the significant correlations in Table 5 demonstrate high consistency in the quality of E.C. vs. SO₄^2–, SAR_or vs Na⁺, and HCO₃^–, and SAR_ad vs Na⁺, and HCO₃^–. The opposite is true for SAR_or and SAR_ad, which exhibit high inconsistency with Mg²⁺ [70].

As shown in Table 6, only four parameters (Sal, Ca²⁺, Cl⁻, and CO₃²⁻) showed no significant change in their mean values [31], a result attributable to the fact that recharge from precipitation has little influence on the magnitude of these four parameters [17].

Figure 5 presents 18 sodium bicarbonate water compositions, that is, four fewer samples than in Figure 6. This difference is attributed to the fact that recharge generates Na⁺ dissolution and reconcentration [24]. The opposite is true for the sodium chloride sulfate composition; in other words, Figure 5 presents 8 samples, three more than Figure 6. This effect is attributed to two causes; a) the precipitation of less soluble ions during the low-flow period [24], and b) the mobilization and redistribution of Cl⁻ of agricultural [71] or marine [72] origin.

The MLR (Equations 6 and 8–10) efficiently predict water quality because, as shown in Table 7, they satisfy the conditions of normality, homoscedasticity (absence of heteroscedasticity), homogeneity, independence (absence of autocorrelation), and minimal multicollinearity, while demonstrating reliable performance [39,42,48,53]. As stated by [39,42], MNR (Equations 7 and 11–14) demonstrate reliable performance.

The Sal models (Equations 15 and 16) showed a good fit, as they met the condition Sr > ?Scr? [42].

The quality models (Equations 6–14) might not be highly efficient, as only a single year of measurement is considered, for both spring and winter. It is recommended to evaluate quality over multiple years to attempt to capture the potential occurrence of extreme variability in precipitation, temperature, and evaporation events in Sinaloa [73]. This extreme variability can destabilize [74] the results regarding observed vs. predicted quality, which are shown in Figure 7 of this study. Another limitation of the models is the complexity and uncertainty inherent in non-linearity [75], as well as the potential failure to select the best-performing function [76].

5. Conclusions

To classify and predict irrigation water quality in spring and winter and based on the relationship of E.C. vs. SAR_or, SAR_co, and SAR_ad, 17 physicochemical parameters were measured in 29 well-water samples. The study was conducted in the coastal portion of the municipality of Guasave, Sinaloa, Mexico.

A non-alkaline pH results in good nutrient availability, which is a major factor in all CO₃²⁻ measurements being zero. The decreases in E.C., TDS, Temp, Salinity, and SO₄²⁻ during winter are due to dilution and a reduction in the depth of s.l. The increases observed in winter for the parameters evaluating sodicity (Na⁺, Ca²⁺, Mg²⁺, HCO₃^–, SAR_or, SAR_co, and SAR_ad) were due to mineral dissolution, leaching of soil salts, cation exchange, and mixing with seawater.

Approximately three-quarters of the total samples showed a high salinity hazard (C3 classification); in other words, the use of these waters is recommended only for irrigating alfalfa, wheat, oats, and barley. To avoid a high risk of sodification (S4 classification), it is suggested to use the water only in winter.

The decline in quality depends significantly on the evaporative effect, minimal recharge, soil leaching, cation exchange, and salt dissolution. High quality consistency is evident when four dependencies are present: 1) E.C. vs. SO₄^2–, 2) SAR_or vs. Na⁺, 3) SAR_co vs. Mg²⁺ and HCO₃⁻, and 4) SAR_ad vs. Na⁺. Mg²⁺ is significant for the occurrence of two major quality inconsistencies related to sodicity (SAR_or and SAR_ad).

The means for Sal, Ca²⁺, Cl^–, and CO₃^2– are not significantly sensitive to dilution from seasonal summer recharge.

91.34% of the total hydrochemical compositions were sodium bicarbonate and sodium chloride sulfate types; spring (18 and 8 compositions, respectively) and winter (22 and 5 compositions, respectively).

Nine predictive models (four MLR and five MNR) of the water quality of wells intended for irrigation are provided. Predicting the rate of change in quality (Equations 12–14) for spring and winter is a methodology that has not previously been applied for the state of Sinaloa. This methodology can encourage the development of early strategies for soil and crop protection against the impacts of irrigation with saline and sodic water.

All models are mathematically validated; however, their predictive capacity could diminish with the occurrence of possible extreme weather events.

This methodology can be applied to any region of the world, and can thus be used to help increase crop yields without compromising the structural stability of the soils.