Integrated modeling of surface water and groundwater systems: a multivariate stochastic perspective

The challenge of using short-term climate data series stems from their limited temporal coverage, which often fails to capture the full spectrum of climatic variability and extremes. This limitation hinders the accurate assessment of long-term trends and the modeling of hydrological processes essential for water resource management and agricultural planning. As a result, there is an urgent need to generate synthetic stochastic climate series that can replicate the statistical properties of observed data over extended periods. This study developed an integrated modeling framework that combines stochastic rainfall generation, monthly evapotranspiration estimation (MASHWIN), and the coupling of the rainfall–runoff model HBV with the groundwater-flow model MODFLOW. The proposed MASHWIN–HBV–MODFLOW approach incorporates uncertainty in precipitation patterns and simulates surface-water and groundwater interactions. The framework was applied to the Aguascalientes Valley aquifer. The stochastic model generated 100 synthetic climate series, which were successfully coupled with the hydrologic and subsurface models. The resulting simulations reproduced the statistical behavior of the observed climate data and enabled analysis of recharge and groundwater-head responses under multiple equiprobable scenarios. The integrated framework provides a comprehensive basis for evaluating hydrological uncertainty and supports long-term groundwater-resource assessment and management in data-limited aquifers.

One of the greatest challenges in hydrogeological studies is estimating recharge, as it cannot be measured directly and must be quantified through indirect methods (Chemingui et al., 2015). The recharge rate varies based on land type and meteorological conditions, leading to significant variability in the values obtained (Águila et al., 2019). Recharge is a fundamental component of hydric balances at the basin level and is essential for groundwater modeling. In semiarid or arid regions, estimating recharge is even more challenging due to limited precipitation, high rates of evapotranspiration, and changes in land use (Hernández-Marín et al., 2018).

The choice of an appropriate technique for quantifying recharge depends on various practical considerations, including spatial and temporal scales, data availability, and the objectives of the study (Chemingui et al., 2015). Mathematically modeling recharge in aquifers also requires precise determination of various components of the hydrological water cycle, such as superficial runoff and evapotranspiration. Numerical models that integrate both superficial and subsurface processes can help address this issue (Ware et al., 2023). However, the process of converting the volume of precipitation into recharge volume remains poorly understood (Kotchoni et al., 2019).

Modeling presents several challenges, the most significant of which is the accessibility of high-quality data that is essential for constructing reliable models. While various hydrological variables have been measured since the 19th century, data collection in developing countries has occurred over a shorter timeframe (Benavides-Bravo et al., 2015; Cárdenas et al., 2017). Information on climatic variables is crucial for calibrating hydrological models (Ghale, 2019; Osuch et al., 2019).

Due to the limited information available in historical databases (such as surface runoff and piezometric data), stochastic models must extend simulation times and propose probabilistic scenarios (Beven, 2021; Nguyen et al., 2023).

Stochastic rainfall patterns have a significant impact on both surface runoff and groundwater recharge, especially in the context of climate change scenarios (Oyati and Yahaya, 2017; Gao et al., 2020; Mishra, 2023). A monthly stochastic weather generator is designed to simulate weather data by producing synthetic time series that reflect the statistical properties of observed climatic variables (Semenov and Brooks, 1999; Ebtehaj et al., 2020; Su et al., 2020). Different authors have proposed multisite multivariate stochastic models for generating synthetic series. These include the Periodic Autoregressive model (Edwin and Martins, 2014; Tsoukalas et al., 2018), data normalization through various transformations (Caloiero et al., 2016; Ebtehaj et al., 2020), the implementation of a Neymann–Scott Rectangular Pulse (NSRP) model (De Luca and Petroselli, 2021), a copula-based conditional quantile model (Wang et al., 2023), and a multilevel multi-objective stochastic approach (Zhang et al., 2020).

Various hydrological models have been used alongside stochastic runoff generators (including the work of Suárez-Almiñana et al., 2017), climate generators paired with hydrologic models (Downs et al., 2023; Jorquera and Pizarro, 2023), stochastic generators combined with the GR4J model (Nguyen et al., 2023), the TUW model (MacDonald et al., 2024), and the WaterGAP3 model (Forte and Rossi, 2024). Additionally, the hyper-multimodel framework (de Assis de Souza Filho et al., 2023) has also been employed. On a different note, several runoff-groundwater simulation models have been developed, such as SWAT-MODFLOW (Kim et al., 2008; Bailey et al., 2016), HBV-MODFLOW (Ejaz et al., 2022), GFLOW-MODFLOW (Feinstein, 2012), Qbox-MODFLOW (Hashemi et al., 2015), and WetSpa-MODFLOW (Tam and Nga, 2018).

Despite the increasing use of stochastic climate generators and coupled surface groundwater models, an important methodological gap remains in the explicit propagation of multivariate climate variability through rainfall runoff and distributed groundwater flow models. Existing integrated frameworks such as SWAT-MODFLOW primarily focus on the spatial exchange of hydrological fluxes between surface and subsurface domains (Bailey et al., 2016; Kim et al., 2008), while HBV-based geohydrological approaches have demonstrated the potential of conceptual rainfall runoff models for long-term groundwater predictions under reduced data and computational requirements (Ejaz et al., 2022). However, these approaches do not explicitly address the propagation of multiple stochastic meteorological realizations from a multivariate weather generator through a rainfall runoff model and into a spatially distributed groundwater-flow model.

This distinction is particularly relevant for semi-arid aquifers, where limited observational records and strong temporal variability in precipitation can introduce substantial uncertainty in recharge estimation. Previous studies in Mexican semi-arid aquifers have addressed groundwater recharge and flow using water-table fluctuation, water-balance, GIS-based, and numerical groundwater approaches, including applications in the Aguascalientes valley (Hernández-Marín et al., 2018). However, these studies have not explicitly propagated multivariate stochastic climate realizations through a rainfall–runoff model and subsequently into a distributed groundwater-flow model. The present study addresses this gap by introducing MASHWIN as an uncertainty-generation layer within the MASHWIN–HBV–MODFLOW framework. The main contribution of this study is the development of a sequential uncertainty-propagation framework that integrates a multivariate stochastic weather generator (MASHWIN), a conceptual rainfall–runoff model (HBV), and a distributed groundwater-flow model (MODFLOW). Unlike previous HBV-based geohydrological approaches and conventional surface–groundwater coupling frameworks, the proposed methodology generates 100 alternative meteorological realizations and propagates their effects sequentially from climate forcing to runoff and recharge and, finally, to groundwater-level responses. This structure provides a quantitative assessment of how climatic variability propagates through the hydrological system and contributes to groundwater uncertainty in data-constrained semi-arid aquifers. To achieve this objective, the study addresses the following research questions: (a) How does the stochastic variability in precipitation and evapotranspiration influence recharge estimates derived from a rainfall–runoff model? (b) To what extent does the variability in recharge affect groundwater level responses? (c) Can stochastic climate series serve as a reliable substitute for observed data in groundwater modeling applications, especially under data-scarce conditions?

The main contribution of this study is the development of a sequential uncertainty-propagation framework that integrates a multivariate stochastic weather generator (MASHWIN), a conceptual rainfall–runoff model (HBV), and a distributed groundwater-flow model (MODFLOW). Unlike previous HBV-based geohydrological approaches and conventional surface–groundwater coupling frameworks, the proposed methodology generates 100 alternative meteorological realizations and propagates their effects sequentially from climate forcing to runoff and recharge and, finally, to groundwater-level responses. This structure provides a quantitative assessment of how climatic variability propagates through the hydrological system and contributes to groundwater uncertainty in data-constrained semi-arid aquifers.

Study area

The Aguascalientes Valley aquifer (Figure 1) is one of the 653 administrative aquifers in Mexico, as designated by the National Water Commission (CONAGUA, 2018). Located in Aguascalientes State, this aquifer supplies approximately 94% of the state’s potable water demand (Hernández-Marín et al., 2018). The region has limited natural surface water sources due to its semiarid climate. The average annual precipitation is 510 mm, while the potential yearly evaporation reaches 2,010 mm, with temperatures between 18 and 22 °C (CONAGUA, 2020). The rainy season begins in May and lasts until October, with the highest rainfall typically occurring in July and August.

Aguascalientes State has a total population of nearly 1.5 million residents, with around 60% living in the state capital, Aguascalientes City (INEGI, 2021). The capital is located within the Aguascalientes Valley and is the center for 95% of all economic activities in the state (Sainz-Santamaria and Martinez-Cruz, 2019).

Geologically, the aquifer is located within the Aguascalientes graben, which is bordered by north–south trending normal faults (CONAGUA, 2020). Given that few reliable dating analyses have been performed, the stratigraphy of the graben remains poorly defined and is primarily established through lithological correlations. However, available drilling logs indicate that the aquifer comprises a sequence of sand and gravel sediments, with only minor intercalations of silt and clay (Pacheco-Martínez et al., 2013). According to the Mexican Geological Service (SGM), the general geology of the region is predominantly composed of alluvial material, with conglomerates (specifically sandstone) and tuff found in the southernmost part of the aquifer (SGM, 2005).

The aquifer system is mainly unconfined, functioning as a water table-type aquifer. However, it exhibits confined aquifer hydraulic behavior in certain areas (Hernández-Marín et al., 2018). The aquifer spans an area of 3,129 km 2. Over the past few decades, water levels have significantly declined due to excessive pumping needed to satisfy the increasing demand for water resources in agricultural, municipal, and industrial sectors (Pacheco-Martínez et al., 2013). The estimated extraction rate from the aquifer is 348 million m 3 per year, whereas the natural recharge rate is only about 250 million m 3 per year. This difference indicates a severe overdraft of nearly 100 million m 3 per year, affecting the aquifer’s nonrenewable storage (CONAGUA, 2020). Extraction is considered the sole outflow from the aquifer, which primarily exhibits a north–south flow pattern (CONAGUA, 2015, 2020).

General methodology

The methodology used in this work is shown in Figure 2, with a more detailed overview in Figure 3. First, meteorological variables, specifically rain and evapotranspiration (ETP) series, were collected and processed. These data were then input into a rainfall-runoff hydrological model, which was calibrated using historical runoff data. Moreover, stochastic series were generated based on the initial rainfall and evapotranspiration (ETP) series. These new data series were input into the calibrated hydrological model to obtain infiltration rates, which were used in the groundwater model. The groundwater model was developed after collecting and processing hydrogeological variables and calibrated using observed head data. Finally, multiple groundwater scenarios were generated through simulations of the groundwater model, using the infiltration rates determined by the rainfall-runoff model with the stochastic rainfall and ETP series.

Meteorological data and stochastic series generation

Rainfall and temperature data are collected from meteorological stations (MSs) located in and around the study area (Figure 1). Additionally, hydrometric data are obtained from hydrometric stations (HSs). The information is sourced from databases managed by the Mexican National Meteorological Service (SMN, 2021) and the National Water Commission (CONAGUA, 2016). The ETP series were calculated utilizing the Hargreaves-Samani method (Hargreaves and Samani, 1985), which relies on maximum and minimum temperature series.

The model for hydrological series analysis (MASHWIN) is a multisite multivariate stochastic model developed at the Technical University of Valencia. It is available at: https://aquatool.webs.upv.es/aqt/otro-software/. This model combines multivariate autoregressive and moving average (ARMA) modeling with spatial disaggregation to provide a detailed understanding of hydraulic system behavior. This study proposes using MASHWIN to generate synthetic time series for precipitation and ETP. The model can create monthly synthetic series with various normalization and periodic parameters, as well as calibrate multivariate parameters. Furthermore, MASHWIN enables the simulation of realistic climatic scenarios that replicate the statistical characteristics of historical series. For precipitation and potential evapotranspiration (ETP), stochastic generation was performed using a periodic first-order multivariate AutoRegressive model, AR(1), at a monthly time step. Seven meteorological stations were considered for each variable. The first-order formulation was selected following the principle of parsimony because, for a seven-dimensional system, each autoregressive lag introduces 7×7 coefficient matrix. Therefore, AR(1) requires 49 lag coefficients for each monthly parameterization, whereas MAR(2) would require 98, substantially increasing model dimensionality relative to the available 29-year observational record. The selected structure therefore provides a compromise between preserving month-to-month persistence and limiting parameter-estimation uncertainty and overfitting. Prior to model fitting, precipitation was subjected to a logarithmic transformation to reduce positive skewness, whereas ETP was modeled without a normalization transformation. No spatial disaggregation was performed because the historical and synthetic series corresponded to the same meteorological stations and the same monthly temporal resolution; consequently, no spatial weighting coefficients were required in this application. Although MASHWIN includes spatial-disaggregation capabilities, this option was not activated in the present study. A 5-year stochastic warm-up period was applied before retaining each generated realization. This period was used only to minimize the influence of the initial conditions of the autoregressive process and was discarded before the synthetic series were transferred to the subsequent hydrological modeling stages. It should therefore not be interpreted as a climate-warming scenario or as the warm-up period of either HBV or MODFLOW. After removal of the stochastic warm-up, 100 equiprobable synthetic realizations of 29 years each were retained for precipitation and ETP, corresponding to 348 monthly values per realization. Model performance was evaluated through the comparison of monthly statistics, autocorrelation functions, and cross-correlation matrices between the stochastic series and the statistical characteristics represented by the historical records.

HBV rainfall-runoff modeling and coupling with MODFLOW

The HBV rainfall-runoff model uses rain and ETP series as inputs. Based on these variables, the model generates a total runoff series, which is used in the calibration process. Additionally, the model produces other outputs, such as infiltration rates, which serve as recharge series for the groundwater model.

The calibration of the HBV model was conducted using a small watershed basin generated via the QGIS GRASS tool (GRASS Development Team, 2017), with HS ID 12239 providing the historical observed runoff series (Figure 1). Because sufficient independent runoff observations were not available for basin-specific calibration of AA and MA, the calibrated HBV parameter set obtained from the gauged calibration basin was transferred unchanged to both subbasins (Hughes-Lomelín, 2023). Differences in simulated recharge between AA and MA therefore reflect differences in basin-specific meteorological forcing under a common HBV parameterization. The infiltration outputs obtained from the calibrated HBV model were subsequently used as recharge inputs for the groundwater model, as described in Section 2.2.4.

The coupling between the rainfall–runoff model (HBV) and the groundwater flow model (MODFLOW) was achieved by transferring infiltration-derived recharge fluxes. In the HBV model, infiltration refers to the portion of precipitation that percolates through the soil and contributes to groundwater recharge. This variable was extracted from the calibrated HBV simulations and used as the primary recharge input for the groundwater model.

The HBV model was applied at the subbasin scale, distinguishing between the Alto Aguascalientes (AA) and Medio Aguascalientes (MA) basins. For each subbasin, a spatially lumped infiltration time series was generated. These recharge values were then assumed to be spatially uniform within each corresponding subbasin and were distributed across the MODFLOW grid cells located within the basin boundaries. This assumption aligns with the regional scale of the study and the limited spatial resolution of available hydrological data.

Temporally, the HBV model operates on a monthly time step, whereas the MODFLOW model was configured using annual stress periods to maintain consistency with the temporal resolution of the available groundwater abstraction, boundary-flow, and piezometric information. Monthly HBV infiltration values were therefore aggregated using a hydrological year extending from October of the preceding year to September of the corresponding year; i.e., October 1985 through to September 1986, with annual recharge assigned to year 1986. Thus, the annual recharge applied to MODFLOW was calculated as the sum of the twelve monthly infiltration values, preserving the total recharge volume associated with each hydrological year while smoothing its intra-annual timing.

This coupling approach allows for the transmission of climatic variability, including stochastic precipitation inputs, through the hydrological system and into the groundwater model. This enables an evaluation of how uncertainty in recharge affects groundwater levels.

MODFLOW model configuration and calibration

The groundwater flow model was developed using MODFLOW within the ModelMuse interface, based on a simplified conceptualization of the Aguascalientes Valley aquifer that aligns with previous studies and the availability of regional data. This model builds upon an earlier framework that was developed and evaluated for the same aquifer system (Hughes Lomelin et al., 2025).

The model domain was discretized into a structured grid consisting of 828 cells (2 × 2 km), arranged in 46 rows and 18 columns, with 335 of these defined as active cells representing the aquifer system.

The aquifer was modeled as a single-layer system with variable thickness, ranging from approximately 300 meters in the western region to 400 and 450 meters in the eastern region, based on available geological information (Pacheco-Martínez et al., 2013; Hernández-Marín et al., 2018). Although the aquifer is known to exhibit heterogeneity and multilayer behavior, a simplified single-layer representation was chosen to effectively capture the predominant regional groundwater flow dynamics, given the scale of analysis and the limited availability of detailed hydrostratigraphic data.

Hydraulic conductivity values were spatially distributed throughout the model domain based on lithological information and interpolation of available data. Storage parameters, including specific yield and storage coefficient, were assigned using values reported in regional studies (CONAGUA, 2020) and were adjusted during the calibration process.

The model was implemented under transient conditions over a 29-year simulation period from 1985 to 2014, using annual stress periods. The initial piezometric surface was defined based on observed groundwater levels corresponding to the starting year of the simulation.

The model considers key inputs and outputs, including vertical recharge, groundwater extraction, and lateral groundwater flow. Vertical recharge was derived from infiltration outputs obtained from the HBV rainfall–runoff model, as described in Section 2.2.4. Groundwater extraction was incorporated using annual pumping rates based on officially reported water concessions https://sigagis.conagua.gob.mx/locrepda20/. The model boundaries were defined based on the hydrogeological limits of the aquifer, and lateral groundwater flow was represented as a constant horizontal flux across these boundaries. This approach is consistent with regional groundwater flow patterns described in previous studies and reflects the exchange of groundwater on a regional scale.

This conceptualization effectively illustrates the primary hydrodynamic processes governing groundwater flow within the aquifer while maintaining a complexity level consistent with the available data. A more in-depth description of the model’s construction, parameterization, and calibration procedure, including calibration and validation against observed piezometric data, can be found in Hughes-Lomelín (2023) and Hughes Lomelin et al. (2025).

The groundwater model was calibrated by comparing the observed piezometric (or head) water levels with the simulated head values generated by the model. Comparisons were made at specific grid cells (Figure 1) selected based on their proximity to observation wells that had the highest density of data points, primarily distributed along the north–south axis of the aquifer. The model was executed and calibrated using recharge inputs obtained from a deterministic superficial rainfall-runoff model, which utilized actual (unsimulated) rain and ETP values.

During the calibration process, several parameters were adjusted, including hydraulic conductivity, storage coefficient, and specific yield. Model performance was assessed through various statistical metrics commonly referenced in the literature (Boyce et al., 2015; Radmanesh et al., 2020; Roohollah et al., 2020; Akter and Ahmed, 2021; Zanini et al., 2021; Mini et al., 2022), as summarized in Table 1. A split-sample strategy was used, where the observed head data were divided into two periods: the first 80% was allocated for calibration, and the remaining 20% was designated for validation. Performance was evaluated separately, with no further parameter adjustments made during the validation phase.

Because the available piezometric observations extend only through 2014, calibration and validation were restricted to the historical 1985–2014 period. Groundwater simulations after 2014 therefore represent stochastic scenario projections rather than independently validated forecasts.

Generation of multiple stochastic groundwater scenarios

The stochastically generated rain and ETP series were processed through the calibrated rainfall-runoff model to produce the corresponding vertical recharge series (ensemble). These recharge series were then input into the calibrated groundwater model, resulting in a total of 100 scenarios. These scenarios were evaluated and compared to assess differences arising from variability in recharge. Throughout the evaluation period and across all scenarios, other inputs, such as groundwater extraction and lateral recharge from adjacent aquifers, were kept constant.

The multisite multivariate stochastic model generated 100 synthetic series, each with a length matching the historical series (29 years). Precipitation was parameterized monthly using logarithmic normalization and included a 5-year warming period. For potential ETP, a normal distribution was assumed within 95% confidence limits.

The generated series are presented in Figures 4, 5, where the monthly mean values for ETP and rainfall can be observed. The synthetic ETP series were adjusted for the various stations analyzed, exhibiting a variability of approximately 10 mm/month compared to the historical series. Figure 4 shows that the months from March to May experience the highest ETP across all stations. In contrast, ETP during winter months averages around 100 mm, increasing in spring and decreasing in summer and autumn.

Precipitation in the study area reflected the highest values in the generated series but was adjusted to align with the average historical period. At the different stations, there is variability of 60 mm/month during the rainy months, which decreases to 10 mm/month in the dry season. The equiprobable series generated are adapted to match the monthly historical series.

Recharge series

One hundred stochastic time series of rainfall and ETP were generated and processed through the calibrated HBV model, resulting in monthly recharge series for the Alto Aguascalientes (AA) and Medio Aguascalientes (MA) basins, which were then aggregated by hydrological year (October–September) and expressed as annual recharge depth series. These are summarized in Figure 6, which also compares the range of simulated recharge values with the deterministic output from the HBV model.

Figure 6 shows that the AA basin experiences lower and more stable recharge compared to the MA basin. The median annual recharge depth in the AA was approximately 2 mm throughout the study period, with relatively narrow interquartile ranges. This indicates limited variability in the hydrological response. In contrast, the MA basin exhibits greater fluctuations in recharge depth, generally ranging between 4 and 6 mm per year, with a wider dispersion in values and noticeable high outliers, particularly during the years 2018 and 2035. While overall recharge is lower in the AA basin, there is a tendency for a greater number of outlier values compared to the MA basin. This suggests that occasional extreme events, whether unusually wet or dry, may result in values falling outside the expected distribution.

The deterministic HBV series primarily aligns closely with the ensemble median for both basins. However, slight deviations during peak years suggest that individual stochastic simulations may occasionally diverge from average ensemble behavior.

To provide physical context for the contrasting recharge responses, the spatial distributions of land use/land cover, lithology, relative soil permeability, and terrain slope were compared between the AA and MA basins (Figure 7).

First, AA covers approximately 2,791 km 2, whereas MA covers approximately 801 km 2. This difference in basin area affects the comparison of recharge expressed as absolute volume and as equivalent water depth. For the historical deterministic HBV series, mean annual infiltration volume was greater in AA (2.36 hm 3 yr. −1) than in MA (1.80 hm 3 yr. −1). However, after normalizing these volumes by basin area, the corresponding mean annual recharge depths were approximately 0.85 mm yr. −1 for AA and 2.25 mm yr. −1 for MA. Thus, although AA generates a greater absolute infiltration volume, recharge per unit area is greater in MA.

Agriculture is the dominant land-cover class in both basins, occupying approximately 53.9% of AA and 49.3% of MA, while built-up areas represent approximately 6.6 and 2.0%, respectively. Conversely, forest, shrubland, grassland, and secondary shrub vegetation collectively occupy approximately 39.0% of AA and 48.2% of MA. Relative soil permeability is predominantly low in both basins, accounting for approximately 66.9% of AA and 73.0% of MA; medium-permeability soils occupy approximately 27% of each basin, whereas high-permeability soils account for approximately 5.1% of AA and are essentially absent from MA. The lithological composition differs more markedly. AA is characterized primarily by alluvium (41.9%), rhyolite–ignimbrite (27.7%), and sandstone–polymictic conglomerate (13.6%), whereas MA is dominated by sandstone–polymictic conglomerate (40.2%), rhyolitic tuff (28.7%), and rhyolite–ignimbrite (21.3%). In contrast, terrain slope is similar between the two basins, with mean values of 2.37° in AA and 2.18° in MA and median values of 1.18° and 1.13°, respectively.

An examination of the original monthly HBV infiltration series confirmed a marked intra-annual concentration of recharge. Across the historical period, approximately 66.3 and 65.8% of total simulated infiltration in the AA and MA basins, respectively, occurred between May and October. Moreover, the highest-infiltration month within each hydrological year accounted for a median of 43.6% of the corresponding annual total in AA and 41.0% in MA. These results indicate that annual aggregation removes substantial seasonal structure from the recharge signal, although the total recharge volume assigned to each annual MODFLOW stress period is conserved.

Groundwater model performance

The calibrated groundwater model, which utilized recharge series obtained from the deterministic HBV model, reached the following average values for the fitting metrics (Table 2).

The RMSE indicates consistent model behavior across observation wells, as reflected in the low standard deviations. The NSE values demonstrate an excellent agreement between the observed and simulated head series, while the R values confirm a strong linear relationship between the modeled and observed values. During the validation phase, the model exhibited slightly lower performance compared to the calibration phase. Nevertheless, the metrics from both processes indicate that the groundwater model is well-calibrated and capable of reliably reproducing head patterns.

Groundwater scenario results under stochastic recharge forcing

The 100 stochastic recharge series obtained from the HBV model were input into the calibrated MODFLOW groundwater model to assess the impact of recharge variability on groundwater levels. Since all other boundary conditions—such as extraction rates and lateral inflows from adjacent aquifers—were kept constant across scenarios, the results reflect only the influence of recharge. The groundwater level time series and spatial patterns were analyzed for temporal trends, variability, and aquifer responses to these fluctuations in recharge.

Figure 8 shows the groundwater levels across the 12 calibration cells.

Over the 29-year simulation period, the cells demonstrate a consistent declining trend, highlighting the imbalance between groundwater extraction and recharge. The deterministic simulation (shown in blue) aligns closely with the median of the ensemble range, indicating that the stochastic recharge series effectively represents a plausible average hydrological condition.

Figure 9 illustrates the spatial and temporal evolution of groundwater levels under stochastic conditions.

A general decline in water levels is observed across most areas of the aquifer during the simulation period, particularly in the northern and central regions, which correspond to the primary extraction zones (agriculture in the north and urban-domestic use in the center). The median response of the ensemble closely matches the deterministic model. However, as shown in Figures 8, 9, the ensemble divergence increases over time. To complement the graphical comparison, the agreement between the deterministic HBV–MODFLOW groundwater-head trajectory and the central MASHWIN–HBV–MODFLOW ensemble response was quantitatively evaluated using RMSE, NSE, and Pearson’s correlation coefficient (R) (Table 2). Across the aquifer evaluated cells, RMSE, with an average of 21.33 m. NSE mean 0.87, while R mean 0.95. These results indicate an extremely high correspondence between the deterministic trajectory and the central tendency of the stochastic simulations. Nevertheless, this agreement should not be interpreted as an absence of uncertainty, because the performance statistics characterize only the central ensemble response. The boxplot dispersion represents the variability among the 100 individual stochastic realizations and shows that recharge uncertainty propagates into an increasingly broad range of possible groundwater-head responses during the simulation period.

Rationale for selecting MASHWIN and methodological trade-offs

The selection of the stochastic generator is an important methodological decision because the characteristics of the synthetic meteorological series directly affect the uncertainty propagated through the HBV and MODFLOW models. In this study, MASHWIN was selected based on the requirements of the complete modeling framework rather than assuming that it is universally superior to other stochastic approaches. MASHWIN provides a multivariate stochastic framework based on autoregressive-moving-average formulations for generating synthetic monthly series while representing relevant statistical and temporal dependence characteristics of the historical record. This is particularly suitable for the present study because the generated meteorological realizations are subsequently used as inputs to HBV and their effects are propagated to groundwater recharge and MODFLOW simulations. Thus, the stochastic component serves as a mechanism for transferring climatic variability into uncertainty in hydrological and groundwater responses.

Several alternative approaches provide important advantages for specific applications. Neyman–Scott Rectangular Pulse (NSRP) models are particularly appropriate for representing rainfall occurrence, storm clustering, duration, intensity, and temporal dependence, making them attractive for event-based or sub-daily rainfall simulation (De Luca and Petroselli, 2021). However, the present study focuses on monthly meteorological variability and its propagation through a rainfall–runoff and groundwater modeling chain rather than on reproducing individual storm events. SPARTA (Stochastic Periodic AutoRegressive To Anything) provides greater flexibility in representing arbitrary marginal distributions while preserving periodic dependence structures (Tsoukalas et al., 2018). Similarly, copula-based models allow flexible representation of nonlinear and potentially tail-dependent relationships among variables (Wang et al., 2023). These characteristics can be advantageous when the main objective is to reproduce complex distributional or dependence structures. However, they require additional specification and validation of marginal distributions and dependence structures, increasing model-selection and parameter-estimation requirements.

The main limitation of this choice is the monthly temporal resolution, which limits the representation of intra-storm dynamics and therefore makes MASHWIN less suitable for applications focused on sub-daily rainfall extremes or flood peaks. In addition, the synthetic realizations remain conditioned on the statistical characteristics and representativeness of the historical record. Therefore, the 100 generated series should be interpreted as plausible stochastic representations of historical climatic variability rather than deterministic predictions of future extreme events.

Within the proposed MASHWIN–HBV–MODFLOW framework, the principal advantage is that multiple meteorological realizations can be propagated through the hydrological and groundwater models, allowing the effects of climatic variability on recharge and groundwater levels to be quantified. Thus, MASHWIN was selected as a methodological compromise that provides an appropriate balance between multivariate temporal representation, computational feasibility, and compatibility with the objectives of the integrated modeling framework.

Impact of stochastic climate variability on recharge

A key objective of this study was to evaluate how stochastic variability in precipitation and ETP affects recharge estimates derived from the HBV model. The results indicate that the stochastic climate series effectively replicate the statistical properties of the observed data while introducing variability that reflects plausible climatic fluctuations. This is demonstrated by the variability observed in the ensemble recharge series, particularly in the Medio Aguascalientes (MA) basin, where higher variability and extreme values were identified.

The contrasting recharge responses of the Alto Aguascalientes (AA) and Medio Aguascalientes (MA) basins cannot be attributed to a single physiographic or hydrogeological characteristic. The spatial comparison in Figure 7 shows substantial differences in land use, relative soil permeability, and lithological composition, whereas terrain slope is very similar between the two basins. Low-permeability soils dominate both basins and are slightly more extensive in MA, while high-permeability soils occur only in AA. Lithologically, AA is characterized by extensive alluvial deposits, whereas MA contains larger proportions of sandstone–conglomerate and rhyolitic volcanic units. AA also shows a greater degree of agricultural and built-up land use, while natural and semi-natural vegetation is more extensive in MA. Taken together, these patterns indicate that the greater recharge depth and stochastic variability simulated for MA cannot be explained simply by greater mapped soil permeability or gentler terrain. Basin size also explains the apparent contrast between absolute infiltration volume and recharge depth: the substantially larger AA basin produces greater total infiltration volume, whereas normalization by contributing area results in greater recharge depth in MA. Recharge depth therefore provides the more appropriate basis for comparing the relative hydrological response of the two basins.

This interpretation must also consider the structure of the HBV implementation. The same calibrated HBV parameter set was transferred to AA and MA because independent runoff observations sufficient for basin-specific calibration were unavailable. Consequently, the different recharge responses generated by HBV primarily arise from differences in the meteorological forcing applied to each basin and the resulting partitioning of precipitation and evapotranspiration under a common model parameterization, rather than from explicit spatial representation of soil, land-use, or lithological heterogeneity. Nevertheless, the mapped basin characteristics provide important physical context for interpreting these responses. Greater urban and agricultural development in AA can modify infiltration, runoff generation, evapotranspiration, and soil-water storage, whereas the contrasting lithological assemblage in MA may influence subsurface transmission through differences in weathering, consolidation, and secondary permeability. These relationships should therefore be regarded as physically plausible basin-scale controls rather than direct causal mechanisms resolved by the lumped HBV model.

From a modeling perspective, the dispersion observed among the stochastic recharge series demonstrates that a single deterministic simulation captures only one possible hydrological trajectory and does not represent the range of responses associated with climatic variability. The ensemble approach therefore provides additional information on the spread and occurrence of plausible recharge conditions, which is particularly relevant in semiarid and data-scarce environments. This interpretation is consistent with previous studies emphasizing the importance of representing climatic uncertainty explicitly in hydrological modeling frameworks (Her et al., 2019; Lawal et al., 2023; Tang et al., 2023).

Overall, the results support the use of stochastic climate generators to characterize recharge uncertainty beyond the average conditions represented by a deterministic simulation.

Propagation of recharge uncertainty into groundwater levels

To evaluate how variability in recharge affects groundwater levels, stochastic recharge series derived from the HBV model were used as inputs to the MODFLOW groundwater model. The results indicate that variability in recharge leads to a measurable variation in simulated groundwater levels, which increases over time. This effect is particularly noticeable in the widening range of the ensemble simulations compared to the deterministic case.

Despite the variability introduced by stochastic forcing, the mean response of the ensemble closely follows the deterministic simulation throughout the analyzed period. This suggests that, on a regional scale, the deterministic recharge input offers a reasonable representation of average groundwater conditions. However, the spread of the ensemble highlights that individual climatic scenarios can result in significantly different groundwater trajectories, particularly over extended simulation periods.

The progressive widening of the groundwater-level ensemble can be interpreted as a cumulative storage response to differences in recharge under a persistent extraction regime. Because groundwater pumping and lateral inflows were held identical among all scenarios, each stochastic realization modifies the balance between recharge and the prescribed groundwater withdrawals. During sequences of relatively dry years, reduced recharge increases the annual deficit between inflow and extraction, causing additional depletion of aquifer storage and progressively lower simulated heads. Because this deficit can persist over consecutive years, differences among stochastic realizations accumulate rather than being reset at the end of each stress period, producing the widening ensemble observed over time.

Conversely, stochastic realizations containing relatively wet years generate greater recharge and temporarily reduce the imbalance between recharge and pumping, thereby moderating drawdown relative to drier realizations. However, isolated high-recharge years do not necessarily compensate for storage losses accumulated during preceding dry periods, particularly in an aquifer already subject to sustained overdraft. The simulations therefore illustrate a compound hydrological risk: climatic variability controls the magnitude and timing of recharge deficits, while continued groundwater extraction determines how those deficits are translated into changes in aquifer storage and groundwater levels. Although pumping uncertainty itself was not simulated, maintaining the same extraction regime across all realizations isolates how different climatic trajectories can alter groundwater outcomes under persistent anthropogenic pressure. This cumulative behavior is consistent with previous studies showing that uncertainty in climatic inputs can propagate and amplify through groundwater systems, particularly in semiarid regions where recharge is limited (Liu and Gupta, 2007; Yang et al., 2016; Beyene et al., 2023; Akhter and Ahanger, 2025; Wu et al., 2025).

These findings suggest that while deterministic approaches may adequately represent average conditions, they might not capture the full range of potential groundwater responses. Incorporating stochastic recharge scenarios provides a more comprehensive framework for assessing uncertainty, which is essential for long-term groundwater management and planning under variable climatic conditions (Goderniaux et al., 2011; Stakhiv, 2011; Pulido-Velazquez et al., 2015). It should be emphasized that this uncertainty represents interannual and long-term groundwater responses, as the annual MODFLOW stress periods do not resolve seasonal recharge–head dynamics.

Adequacy of the simplified groundwater model

An important question in this study is whether the simplified single-layer groundwater model is sufficient to capture the key dynamics of the aquifer under stochastic recharge conditions. The model demonstrated strong performance during both calibration and validation, as indicated by the statistical metrics (e.g., NSE, RMSE, and correlation coefficients; see Table 2). This suggests that the model can reproduce observed groundwater level patterns with acceptable accuracy at a regional scale.

The model’s ability to reproduce the regional spatial distribution and temporal evolution of groundwater levels indicates that the principal regional head patterns are represented adequately for the objectives of this study. In particular, the model effectively captures the general north–south flow pattern and the observed declining trends associated with sustained groundwater extraction.

The single-layer representation inherently limits the model’s ability to resolve vertical hydraulic gradients, interlayer leakage, confined–unconfined transitions, and differences in storage among hydrostratigraphic units. Although hydraulic conductivity varies spatially in a horizontal manner within the modeled layer, the vertical organization of these properties is represented through effective regional parameters rather than through explicit individual aquifer and aquitard units. This structural simplification can influence how recharge uncertainty is propagated into simulated groundwater heads in more than one direction.

On the one hand, the absence of explicit low-permeability interbeds and vertical hydraulic resistance may cause recharge perturbations to be transmitted more directly to the modeled groundwater head than they would be in a multilayer system. In this sense, the single-layer conceptualization could amplify or accelerate the apparent response to recharge variability by omitting delays associated with vertical percolation and interlayer leakage. On the other hand, representing vertically heterogeneous units through an equivalent regional layer may dampen local responses because contrasting storage properties, vertical head gradients, perched or semiconfined behavior, and localized recharge pathways are averaged within a single simulated head. The direction and magnitude of this structural effect are therefore spatially dependent and cannot be uniquely quantified from the available regional data.

A two-layer sensitivity model was not implemented because the available dataset does not provide sufficient information to independently constrain layer geometry, layer-specific horizontal and vertical hydraulic conductivities, storage properties, interlayer leakance, pumping distribution by depth, or observation-well screen intervals throughout the model domain. Dividing the aquifer into two numerical layers without these constraints would require introducing additional assumed parameters and could create an appearance of greater physical detail without a corresponding increase in model reliability. Accordingly, the good calibration and validation performance of the present model supports its use for reproducing regional groundwater-level patterns, but it should not be interpreted as evidence that vertical flow processes are uniquely represented.

Within these limitations, the single-layer conceptualization is consistent with the objective of this study: evaluating the propagation of stochastic recharge uncertainty through regional and interannual groundwater-level trajectories rather than resolving local vertical flow. The results should therefore be interpreted at this regional scale. Future work supported by detailed hydrostratigraphic information and vertically resolved groundwater observations should test alternative multilayer conceptualizations and quantify the sensitivity of recharge-driven head uncertainty to vertical hydraulic structure.

Transferability of the integrated modeling framework

The applicability of the proposed MASHWIN–HBV–MODFLOW framework to other semiarid and data-scarce aquifers lies primarily in the integration of stochastic climatic forcing with surface-water and groundwater modeling rather than in the direct transfer of the model parameterization developed for the Aguascalientes Valley. Previous applications in Mexico demonstrate the use of hydrological and groundwater models to evaluate climate-driven changes in water availability and spatial and temporal variability in groundwater recharge (Ruiz, 2015; Molina-Navarro et al., 2016). More broadly, data availability and hydrogeological characterization remain recurrent constraints in numerical groundwater modeling across the Americas (Lozano Hernández et al., 2024). Within this context, the contribution of the present framework lies in the sequential integration of stochastic climatic forcing, rainfall–runoff modeling, and groundwater-flow simulation, allowing climatic uncertainty to be propagated through the coupled hydrological system.

Adapting the framework to another basin would require: (1) compilation and quality control of locally available meteorological and hydrological observations; (2) calibration and evaluation of the stochastic climate generator against the statistical characteristics of the observed series; (3) calibration of the rainfall–runoff model using locally available discharge information; (4) definition of spatial and temporal procedures for transferring infiltration-derived recharge to the groundwater model; (5) construction and calibration of a site-specific groundwater conceptual model; (6) propagation of multiple stochastic climatic realizations through the coupled models; and (7) evaluation of the resulting ensemble uncertainty according to the hydrological or management objective of the application. Consequently, the modeling sequence is transferable, whereas MASHWIN and HBV parameters, hydrogeological properties, spatial discretization, stress-period duration, and groundwater boundary conditions must be independently defined according to local characteristics and data availability.

Limitations and future work

The current study has various limitations. First, the recharge estimates depend on assumptions made for the HBV model and the stochastic generation of meteorological variables. As a result, these estimations may not accurately reflect real recharge rates or fully capture long-term climatic shifts or rare extremes. Furthermore, because a common calibrated HBV parameter set was applied to AA and MA, spatial differences in soil, land use, and lithology were not explicitly represented in basin-specific model parameters; future applications with additional runoff observations could incorporate spatially differentiated calibration to better resolve these controls. In addition, the groundwater model assumes static boundary conditions for pumping and land use throughout the simulation period. Although this assumption helps isolate climate-driven variability, it fails to account for trends in urbanization and agriculture, as well as climate change scenarios that significantly impact future water demand and recharge dynamics. Climate and pumping effects may also interact. Increased abstraction occurring simultaneously with sequences of below-average stochastic recharge would be expected to exacerbate groundwater-level decline, whereas pumping reductions could partially offset the effect of dry climatic sequences. Land-use change could further affect this relationship through modifications in infiltration, runoff, recharge, and groundwater demand. Future applications of the framework should therefore combine the stochastic recharge ensemble with multiple pumping and land-use trajectories. Such factorial scenario analysis would allow climatic and anthropogenic uncertainties to be evaluated jointly and would provide a stronger basis for identifying groundwater-management strategies that remain robust across a range of plausible recharge conditions.

A further limitation concerns the temporal extent of the groundwater observations used for model evaluation. The available piezometric dataset extends through 2014, and therefore the post-2014 stochastic simulations cannot be independently validated against observed groundwater heads. This limitation also persists in the most recent official assessment of the Valle de Aguascalientes aquifer published by CONAGUA (2024), which continues to base its piezometric analysis and groundwater-level evolution on the 2000–2014 period. Nevertheless, the updated official water balance continues to indicate substantial groundwater overdraft, providing qualitative support for the persistence of the depletion pressure represented in the simulations. Accordingly, post-2014 model outputs should be interpreted as plausible stochastic groundwater trajectories conditioned on the calibrated historical system and prescribed extraction regime, rather than as deterministic predictions of observed future heads.

The temporal aggregation of monthly HBV recharge into annual MODFLOW stress periods represents an additional source of temporal smoothing in the coupled framework. The monthly HBV results show that infiltration is strongly seasonal, with approximately two-thirds of the simulated infiltration occurring between May and October and individual peak months accounting for a substantial fraction of the annual total. Consequently, the annual groundwater model cannot reproduce seasonal groundwater-level fluctuations, short-term recharge–head lags, or transient responses to individual recharge pulses. Nevertheless, because the monthly values were summed, rather than averaged, the total recharge volume for each hydrological year was preserved. Annual stress periods were also consistent with the temporal resolution of the available pumping, lateral-boundary, and piezometric information. Introducing monthly recharge into the groundwater model without equivalent monthly information for these other stresses and calibration targets would introduce unsupported temporal assumptions. Accordingly, the MODFLOW simulations presented here should be interpreted as evaluations of interannual and long-term groundwater-head trajectories and uncertainty propagation, rather than simulations of seasonal groundwater dynamics.

Future research should delve into integrating dynamic pumping scenarios and land use changes, along with incorporating climate change projections, which are beyond the scope of this work. Improving the spatial resolution of the groundwater model and developing a multilayer hydrostratigraphic representation, supported by layer-specific hydraulic properties, well-screen information, and vertically resolved piezometric observations, would allow future studies to quantify the sensitivity of recharge uncertainty propagation to vertical groundwater-flow processes. Future developments should also evaluate groundwater sensitivity to seasonal recharge variability through shorter MODFLOW stress periods, contingent upon the availability of temporally consistent monthly pumping, boundary-flow, and piezometric records required for model calibration and validation. Future model updates should also incorporate newly available post-2014 groundwater-level observations as they become sufficiently complete and spatially consistent, allowing independent evaluation of the long-term stochastic trajectories generated by the framework.

This study developed and evaluated an integrated modeling framework (MASHWIN–HBV–MODFLOW) to assess how stochastic climatic variability propagates through groundwater recharge and piezometric responses in a semiarid aquifer. The results demonstrate that stochastic climate series provide a useful approach for extending the range of plausible climatic variability represented by limited observational records and for evaluating its propagation through coupled hydrological models.

Stochastic variability in precipitation and ETP produced measurable differences in recharge and groundwater-level trajectories. Although the central tendency of the stochastic ensemble remained close to the deterministic simulation, the dispersion among groundwater scenarios increased through time as recharge differences accumulated under a persistent extraction regime. This demonstrates that a deterministic simulation can represent average behavior while masking hydrological trajectories associated with comparatively dry or wet climatic realizations.

The findings also indicate that the simplified groundwater model adequately reproduces the principal regional head patterns over the historical calibration and validation period and provides a suitable basis for evaluating recharge-driven uncertainty at the regional scale considered in this study. However, the results should not be interpreted as resolving local vertical-flow processes or seasonal groundwater dynamics, which would require finer hydrostratigraphic and temporal information.

For regional groundwater management, the 100 stochastic realizations provide a probabilistic baseline that can complement the deterministic assessment. Rather than relying only on the expected groundwater trajectory, water managers can use the ensemble distribution to identify the probability of comparatively low-recharge conditions, evaluate groundwater-head responses under dry, median, and wet climatic realizations, and identify periods in which climatic variability increases the risk of accelerated groundwater-level decline. In practical terms, the ensemble can support the definition of drought contingency triggers based on the probability of groundwater heads crossing management-defined thresholds, provided that those thresholds are established by the responsible water authority. The framework can also be used to develop probabilistic groundwater-abstraction strategies.

Future model development should therefore follow a prioritized sequence. First, stochastic recharge ensembles should be coupled with dynamic groundwater-extraction scenarios, because pumping is a dominant control on groundwater decline in the Aguascalientes Valley and this extension would provide the most direct management benefit. Second, seasonal groundwater stresses and observations should be incorporated to improve the representation of drought response and short-term management interventions. Third, multilayer hydrostratigraphy, land-use change, and nonstationary climate-change projections should be incorporated to refine spatial and long-term risk assessments. This progression would allow the MASHWIN–HBV–MODFLOW framework to evolve from a climate-uncertainty propagation tool toward a probabilistic decision-support framework for regional aquifer management.