Abstract
Groundwater resources are influenced by a complex interplay of natural processes and anthropogenic activities. In the context of accelerating climate change, the availability of tools able to allow the implementation of robust and evidence-based resource management policies are increasingly necessary. A particularly relevant application of this approach lies in assessing how shifts in agricultural policy may influence water availability. This study presents an advanced unsteady-state numerical groundwater model for Milan metropolitan area, located within the Adda-Ticino basin in Northern Italy. The model, built with MODFLOW-NWT and subjected to a global sensitivity analysis, was calibrated over a five-year period (2014–2018) and was then used to assess the potential effect of irrigation practice changes. Results highlight the substantial role of irrigation in sustaining groundwater levels, with irrigation accounting 57% of average yearly recharge. Scenario simulations show that reductions in irrigation recharge (15, 30, and 50%) could lead to significant groundwater level declines, particularly in agricultural zones surrounding Milan but also extending to the city area. The model is a useful and robust decision-support tool for water authorities, allowing them to assess policy impacts and support the development of sustainable water management plans under changing climatic and anthropogenic drivers.
Highlights
Irrigation recharge from traditional techniques (flooding the fields) plays an important role in the Adda-Ticino basin.
A robust groundwater flow numerical model was developed, allowing to quantitatively assess the groundwater budget in the Adda-Ticino basin.
Scenarios of irrigation practices change lead to significant groundwater depletion, endangering its availability and dependent ecosystems.
1 Introduction
Groundwater depletion is a widespread phenomenon across the globe, especially in recent years (Jasechko et al., 2024). Groundwater management includes the operational decisions and practices that guide abstraction of groundwater, as well as other activities that have an influence on groundwater and the aquifers through which it flows (United Nations, 2022). It is thus a crucial part in avoiding overconsumption while responding to the needs of local communities. Good groundwater governance should enable the conditions for groundwater management to determine who gets which water, how, where and when, leading to implications for society, the environment and the economy. Making decisions in this sector is both important and challenging, under the responsibility of decision makers and administrations.
Good practice in policymaking should be supported by technical and scientific knowledge and tools. In the groundwater sector, the necessary information traditionally comes from the deployment of monitoring wells, with periodical groundwater level measurements and samplings to monitor water quality and quantity. Since the 1960s, with the development of numerical groundwater models (Freeze and Witherspoon, 1966; Zhou and Li, 2011), and lately with data-driven models (Sarma and Singh, 2022), the pool of available tools was enriched allowing to estimate aquifers levels and contaminants concentration in points where no monitoring wells are present, along with the possibility of generating forecasts of future conditions. They provide the opportunity to test and assess the effect of possible policies and management practices on a “model,” a representation of reality, prior to the application in the real world, being able to consider many different aspects and actors involved or possibly impacted. Furthermore, numerical groundwater models are often used to compute the water budget, which quantifies the main inflow, outflow, and storage components of the simulated system and allows the relative contribution of the different flow components to be assessed, supporting groundwater resource evaluation and management. The literature is rich in examples of numerical groundwater models used as decision support and management tools. Numerical models have been successfully used to analyze problems related to regional-scale groundwater availability and sustainability, both by government and academic authors (Zhou and Li, 2011; Islam et al., 2017; Delfini et al., 2024; Alberti et al., 2025; Delsman et al., 2023; Jackson et al., 2011; Feinstein et al., 2010; Brakefield et al., 2015; Piggot, 2002). Their application is very common in contamination management (Antelmi et al., 2020; Xu et al., 2022; Yang et al., 2021; Antelmi et al., 2021b; Chen et al., 2016), as well as in the assessment of saltwater intrusion (Dyring et al., 2023; Coulon et al., 2022, 2024; Alberti et al., 2022; Babu et al., 2018; Izuka et al., 2021; Lotti et al., 2021), application of MAR techniques (Colombo et al., 2024; Dahlke et al., 2018; Perdikaki et al., 2022; Ringleb et al., 2016; Russo et al., 2015), geothermal purposes (Bulté et al., 2021; Omagbon et al., 2021; Antelmi et al., 2021a) or interaction with man-made structures (Sartirana et al., 2022; Zhang et al., 2024).
Various studies on the quantitative and qualitative management of groundwater resources of Milan are available in the literature (Sartirana et al., 2024), but only a limited number have used numerical modeling as an analysis and forecasting tool. Between the late 1990s and early 2000s, several groundwater models were developed for the city of Milan (Beretta et al., 2004; Giudici et al., 2000; Alberti et al., 2000), primarily aimed at investigating the causes of the observed rise in groundwater levels and supporting the design of potential mitigation strategies. These early models were steady-state, two-dimensional structures with coarse spatial discretization, typically limited to the local or municipal scale. It was not until 2014 that a basin-scale model was implemented (Alberti et al., 2016), as part of the PLUMES project coordinated by ARPA Lombardia and Regione Lombardia. Developed using MODFLOW-2005, this model encompassed the entire Adda-Ticino-Po hydrogeological basin, covering the entire Province of Milan. This configuration enabled the calculation of the water balance for Milan’s three principal aquifers (A, B, and C, later explained in section 2.1). From 2016, a growing number of modeling applications for the city of Milan is available in the literature, both in the field of groundwater/underground infrastructures interactions (Sartirana et al., 2022; Colombo et al., 2018; De Caro et al., 2020; Gattinoni and Scesi, 2017) and in the simulation of contaminant plumes (Alberti et al., 2018; Pollicino et al., 2021a; Colombo et al., 2020) and thermal plumes (Previati and Crosta, 2021; Barbieri et al., 2025).
However, no further work on the quantitative aspects of water resources has been presented until 2020. In De Caro et al. (2020), it is considered the same model domain of Alberti et al. (2016), but focusing only on the first two aquifers (shallow unconfined and first semi-confined aquifer, called A and B, see section 2.1) by means of a finite elements model through FEFLOW (Trefry and Muffels, 2007) However, the water budget shows larger fluxes of main components compared to previous models (Alberti et al., 2016), particularly in terms of horizontal groundwater inflows and outflows across the model boundaries. These differences are further confirmed in the most recent study by Alberti et al. (2025), which presents a refined steady-state model of Milan and its surrounding province. The need to improve the understanding and quantification of groundwater resources in the Milan Province requires a continuous update of the research in this field to achieve increasingly reliable tools that are shared across the scientific community. The present paper fits into this line of research, aiming, starting from the work of Alberti et al. (2025), to advance the model’s reliability and evaluate the importance of the groundwater recharge in the context of groundwater mass balance assessment.
Since groundwater is part of the water cycle, surface waters and human usage have to be included in the picture when considering its management (Abbott et al., 2019). European directives, such as the Water Framework Directive (2000/60/EC), encourage management policies that allow water saving since the early 2000s. In the Public Institutions it is widely believed that in any geographical area flooding and furrows irrigation methods are just obsolete and wasting water systems. The European Union has promoted several regulations to improve irrigation efficiency, such as the Regulation (EU) 2020/7412020 on water reuse. More recently, the Common Agricultural Policy has reinforced efforts to enhance water efficiency and reduce water losses (European Commission, 2026). However, these policy decisions could impact groundwater resources: excess irrigation water that does not reach the crops is not wasted, but percolates becoming a resource for aquifers’ recharge and underground water storage. Irrigation through these methods represent a considerable portion of groundwater recharge.
Surface irrigation systems typically exhibit field application efficiencies ranging from approximately 40 to 65% (Ali et al., 2025; Howell, 2003), indicating that about 35 to 60% of the applied irrigation water infiltrates below the root zone. These values are consistent with the efficiency ranges reported by the FAO for surface irrigation systems (Brouwer et al., 1989). In Northern Italy, field observations by Masseroni et al. (2022) reported water use efficiencies ranging from 16 to 56% under traditional border irrigation practices, due to the high drainage capacity of the soils. Considering the amount of water used by agriculture in agriculturally intensive areas (e.g., 4′600 m3/ha in Italy) (Istat, 2014), irrigation through these methods represents a significant part of groundwater recharge. Giuliano (1995) estimated that the irrigation accounts for about half of the total recharge in the Po basin (more than 7 billion m3/year). Different papers account for the large contribution of the irrigation return flow in the Po plain area, with values ranging between 55–88% of the total recharge, depending on the considered area, period and crop (Rotiroti et al., 2023; Lasagna et al., 2020; Vassena et al., 2012; Baker et al., 2022; Rotiroti et al., 2019; Redaelli et al., 2026).
The development of a tool able to give outcomes of different management choices has thus become necessary when dealing with complex systems. In this work an advanced unsteady-state groundwater numerical model of Milan and its Functional Urban Area [FUA, from OECD (2012)], set in the Adda-Ticino basin, is presented. It aims to support public authorities and water managers to evaluate different alternatives in the decision-making process when they deal with groundwater management. The model is the first advanced numerical tool available for this area with this level of detail, able to successfully represent the recharge dynamics driven by irrigation and water extraction practices. After the calibration and validation processes, the developed model was used to simulate three scenarios of decreasing irrigation recharge, simulating conditions which can mirror blue water scarcity and changes in irrigation practices, quantitatively estimating the variations in groundwater levels. The scenarios are focused on the long period and not on the transitory phase, to get insight of what might happen in the future if the whole irrigation system is switched to more efficient techniques. This allows a comprehensive understanding of the linkages between the irrigation system, groundwater bodies and water availability in the area.
2 Materials and methods
2.1 Study area
The Adda-Ticino basin is located in the Po plain in Northern Italy, characterized by a continental temperate climate, with warm summers and cool and humid winters. The precipitation ranges from 700 to 1′200 mm/yr. The Adda and Ticino rivers are located at the eastern and western boundaries of the study area, respectively, as shown in Figure 1a. The main city in the area is Milan, which with its hinterland constitutes one of the most urbanized and industrialized areas in Europe, with 3.2 million people and a population density of 2040 inhabitants/km2 (Istat, 2023). Nevertheless, Milan and surrounding provinces (Varese, Lodi, Pavia, Monza e Brianza) still present a strong agricultural activity supported by a dense network of irrigation canals. The main branches of these canals, such as Navigli Lombardi (Grande, Bereguado, Pavese and Martesana), Muzza and Villoresi canals (Figure 1a), bring water from the Adda and Ticino rivers to farmlands thanks to a network of minor canals. The crop fields are usually irrigated from late April to September (spring–summer period) mainly through the flooding irrigation technique (Zucaro and Corapi, 2009), while the minor canals are kept dry during fall and winter seasons. The area is also characterized by the presence of the Olona, Lambro and Seveso streams, that flow from the northern boundary of the area towards south, crossing the city of Milan.
Figure 1
The geology of the area is characterized by quaternary alluvial sediments composing a sequence of 4 hydrogeological units of interest beneath the area, called A, B C and D [as defined by Regione Lombardia and AGIP (Regione Lombardia and ENI-AGIP, 2002) and adopted by many authors in literature (Alberti et al., 2025; Sartirana et al., 2024; Gorla et al., 2016; Colombo et al., 2019, 2020; Pollicino et al., 2021b; Bonomi et al., 2010; Gorla, 2001)]. Aquifers A and B are the most exploited both for drinking and industrial purposes, being the most transmissive, and therefore are the only two considered in the present study. Moreover, Aquifer B and C are well separated by a clay layer, that limits the exchanges between the two aquifers and therefore it can be considered an impervious lower boundary. A representative cross-section is shown is Figure 1b while more detailed geological cross-sections are found in Alberti et al. (2025).
Groundwater is extensively exploited in the area. Drinking water is completely supplied by groundwater, that is mainly pumped from Aquifer B, while both aquifers are used for industrial purposes and other private uses. Even if groundwater levels rose since 1970, when many industries ceased their activities and stopped pumping (Alberti and Francani, 2001), management of groundwater is needed to regulate groundwater exploitation and control its quality, as required by Europe through the Water Framework Directive (2000/60/EC) and the Water Safety Plans (WSP, World Health Organization). Within the study area more than 3′000 pumping wells and many surface water bodies (streams, irrigation canals) are present, interacting with the aquifer and influencing its water level over time. Different studies have shown how irrigation recharge in the agricultural areas surrounding the city of Milan plays an important role in supplying the underground water system (Alberti et al., 2016; Alberti and Francani, 2001; Colombo et al., 2024; Giuliano, 1995), but a groundwater numerical model able to represents the groundwater behavior in unsteady-state conditions is still missing.
2.2 Transient numerical model
The numerical model was built using MODFLOW-2005 (Harbaugh, 2005) adopting its Newton–Raphson formulation (MODFLOW-NWT) (Niswonger et al., 2011), due to its robustness and the capability to manage dry cells rewetting. The model domain covers an area of about 1880 km2, including the entire Metropolitan city of Milan and some municipalities of the Monza-Brianza, Varese, Pavia and Lodi provinces. The Western and Eastern boundaries of the domain are represented by the Ticino and Adda rivers, respectively, while the northern and southern boundaries were selected to be far enough from the area of interest (the city of Milan) (Alberti et al., 2025).
The model was built starting from a steady state model elaborated by Alberti et al. (2025), where the conceptual model of the area was built considering about 40′000 stratigraphies, which were used to build a 3D geological model that allowed to identify the different hydrostratigraphic units and subunits (described in section 2.1.1). The results of the steady state model presented in Alberti et al. (2025) were used as initial head condition for the current model. Such model was then set up to simulate five years of system behavior (January 2014–December 2018) for which piezometric measurements were available from ARPA Lombardia, MM S.p.a and CAP Holding S.p.a. The elements capable of influencing the hydrogeological balance (withdrawals, meteoric recharge, irrigation recharge and internal conditions) were used to define the time discretization as discussed later in the manuscript, accordingly to the stresses induced by these elements and are shown in Figure 2.
Figure 2
2.2.1 Spatial and temporal discretization
The model domain was horizontally discretized adopting a uniform grid with square cells of 100 m sides, considered sufficient to ensure good precision and good computational efficiency at the same time, considering the model extent. The cells to the West of the Ticino River and to the East of the Adda River were then inactivated. Vertical discretization was developed to best represent Aquifer A and Aquifer B. Relying on the conceptual model, 7 layers were implemented (Figure 1b). The first layer represents the unconfined aquifer (Aquifer A), the second layer reproduces the hydraulic separation lens (aquitard) present between the two aquifers. Layers from 3 to 7 represent the aquifer sub-units and the clay-silty levels that constitute Aquifer B. The model domain is composed of 1′315’265 active cells (338 rows and 663 columns). Figure 1 shows horizontal and vertical discretization.
The temporal discretization was chosen based on the seasonality of irrigation practices. Each year marks four different Stress Periods (SP). Table 1 shows the length of each SP and the percentage of the average channel flow rates applied in the irrigation period within each district. Over the simulation period of 2014–2018, the model run is characterized by 20 SP (4 for each year) for a total duration of 1′825 days.
Table 1
| SP | Time interval | Days | % of the discharge in the considered SP* |
|---|---|---|---|
| 1 | 1 January– 31 March | 90 | 0** |
| 2 | 1 April – 15 June | 76 | 62 |
| 3 | 16 June – 15 September | 92 | 88 |
| 4 | 16 September– 31 December | 107 | 0** |
Temporal discretization and percentage of channel derivation over the maximum available amount (SP: stress period).
*The original data is the discharge provided to the different canals in an average irrigation period, namely from the 30th of April to the 5th of September. Since the SP discretization splits the period and the SPs boundaries do not exactly match, the discharge has been multiplied by the fraction of actual days vs the modeled days (47/76 and 81/92, respectively). **Navigli and Muzza canals are also active in autumn and winter, so their respective losses have been included.
2.2.2 Boundary conditions
Adda and Ticino rivers represent physical eastern and western limits of the model and are implemented using a Dirichlet (fixed or Constant Head, CH) condition in the first layer. The two rivers show a relatively constant regime because the flow rates are regulated through dams (Miorina for Ticino and Olginate for Adda). The most consistent variations in hydrometric levels are linked to temporally limited episodes. Furthermore, the official piezometric maps elaborated by Città Metropolitana di Milano and Regione Lombardia [Città Metropolitana di Milano (Dataset della categoria AMBIENTE), 2024; Regione Lombardia (Geoportale della Lombardia), 2014] do not show a significant change of the groundwater-river relation over time. For these reasons, it was assumed to keep the river heads constant, therefore Dirichlet conditions were assigned.
A General Head Boundary (GHB, a Cauchy-type or head-dependent flow condition) was implemented for the northern and southern boundaries of the model, representing hydraulic boundaries. In this case there were no physical limits capable of controlling the groundwater level (as in the case of Adda and Ticino); however, the distance of the North and South limits from the area of most interest (city of Milan) was considered large enough to prevent large effects on the model results. GHBs were only defined where the observed piezometric head was higher than the bottom of the layers. The GHB head values were kept constant over time using the head values of 2014 already defined in Alberti et al. (2025) due to limited data availability at the borders of the model and a conductance in the order of 10−3 m2/s (variable from cell to cell) was implemented, according to the aquifer characteristics.
2.2.3 Internal conditions
2.2.3.1 Streams and canals
Aside from Adda and Ticino rivers, the existing streams within the model domain (Olona, Lambro and Seveso) and the main artificial canals used for irrigation purposes (Villoresi canal, Navigli Lombardi and Muzza canal) were represented in the model. The rivers were represented through a Cauchy boundary condition using the MODFLOW RIV package (considering conductance values for the different streams calculated on river geometry and hydraulic conductivity of the riverbed sediments), while the main artificial canals were simulated as water losses by assigning a Neumann (prescribed flow) boundary condition through the MODFLOW WEL package. To evaluate the variations of river regimes in the period considered, the hydrometric levels were analyzed in correspondence of the hydrometers located along their course, within the modeling domain, without observing significant changes in the water levels and in the flow rates, as visible in Supplementary Figure S1. Therefore, the heads inserted through the RIVER package (Cauchy condition) were kept constant over time.
The different irrigation canals that contribute to feeding the groundwater were instead included as estimated distributed losses along their path (considering the canal bed type and the data retrieved from the irrigation district managers). Since the irrigation network is not always active throughout the year, due to the seasonal need of irrigation water, the prescribed flow rate was varied considering the average flow rate flowing through the canals in the various SPs, in compliance with what is reported in Table 1.
2.2.3.2 Groundwater withdrawals
The pumping wells present in the area (3′037) were represented using a Neumann condition by means of the WEL package. The pumping rates varied for each SP from 2014 to 2018, according to data availability. Wells present within the domain are derived from the comparison of different databases:
SIF (Sistema Informativo Falda – Groundwater Information System): Milan’s province database (now dismissed);
CUI (Catasto Utenze Idriche – Inventory of Water Utilities): database of Regione Lombardia;
Water manager databases (Metropolitana Milanese s.p.a., CAP Holding s.p.a., BrianzAcque s.r.l.).
2.2.3.3 Lowland springs (fontanili)
The area is characterized by the presence of lowland springs (locally known as “fontanili”), that were represented through the DRN package (Drain), fixing the stage equal to the bottom of the spring heads (retrieved from a 5×5 m DTM provided by Regione Lombardia) and the conductance equal to 0.2 m2/s, according to the average spring characteristics and sediment hydraulic conductivity.
2.2.4 Recharge
The groundwater recharge was supplied as input to the model as the flow rate of water per unit of area that infiltrates the subsoil and feeds the water table. The value of this flow is given by the sum of the infiltrations due to meteoric precipitation and anthropic activities (irrigation practices and losses of the water networks). The recharge is implemented, using the MODFLOW RCH package, only in the first active layer of the model.
The recharge was estimated through Soil Water Balance 2 (SWB2) (Westenbroek et al., 2018). In the configuration used, the software provides an estimate of the distributed meteoric recharge on a 100×100 m resolution grid, considering the soil water balance in terms of interception, precipitation, evapotranspiration and soil moisture as well as land cover and impervious surfaces. Daily precipitation and temperature (minimum and maximum) were retrieved from the E-OBS database (version 24.0e, Cornes et al., 2018). This data is available at a 0.1° latitude and longitude resolution and was clipped at the model area and in the simulation timeframe (2014–2018). The resolution is deemed sufficient to represent an 1882 km2 area, and the frequent updates to E-OBS permit future utilizations for model updates. The Hargreaves–Samani method was used to estimate the evapotranspiration, as implemented in SWB2 (Hargreaves and Samani, 1985).
The land cover classification was obtained from the DUSAF 2018 database (Regione Lombardia, 2018), which is based on the Corine Land Cover project (Bossard et al., 2000). The model area was thus subdivided into five major classes, mapped in Figure 3: urban, agricultural, green, woodland and pasture, and water bodies. Inside these classes, more accurate land uses from the DUSAF classification were used. To each land use, a curve number was associated to estimate the run-off, following the Curve Number method from USDA Soil Conservation Service (Woodward et al., 2003). The full classification and the associated curve numbers (defined in the SWB input files and corresponding to the different land uses) are provided in Supporting Information (Supplementary Table S1). The impervious surface percentages were selected through literature analysis (Becciu and Paoletti, 2010) and assigned to the urban land use classes, as indicated in Supplementary Table S1.
Figure 3
The irrigation recharge was not estimated through SWB2. In agricultural areas, an estimate of the recharge contribution due to irrigation practice was thus added to the meteoric recharge as a fraction of the canals’ flow rates. Such contribution varies during the year and was computed as a function of the diverted water from the main irrigation canals (Table 1): it was set to zero during fall and winter, while it increased in the spring period and reached its maximum value in the summer period. The values of irrigation recharge were differentiated within the model domain considering the irrigation districts fed by each main canal. For the Villoresi canal, a further distinction has been made considering the different irrigation sub-districts fed by the secondary canals (Figure 3b). The irrigation recharge was calculated by dividing the flow rate distributed through the irrigation canals by the area of the corresponding irrigated district or sub-district, assuming a uniform spatial distribution. The resulting irrigation input was then reduced by accounting for crop evapotranspiration and irrigation efficiency, with the remaining fraction considered as groundwater recharge. The evapotranspiration and plant use contributions were set equal to 30%, a value based on the surface irrigation method’s efficiency, the most common method in the area considered, retrieved from literature (De Pascale et al., 2011) and an additional 5% surface runoff losses. An additional contribution was added in urban areas to account for aqueducts and sewerage losses, calculated for each individual municipality and distributed uniformly over urbanized areas. This percentage was set equal to 12.5% of the groundwater extractions by the water managers in all municipalities, being groundwater the drinking water source in the study area.
2.2.5 Hydrogeological properties
The initial hydraulic conductivity cell-by-cell distribution used in the transient model was retrieved by the work of Alberti et al. (2025), where a steady state model for the same area was previously developed. Then, the hydraulic conductivity distribution was adapted through a history-matching process using the optimization code PEST_HP (Doherty, 2020) using the Pilot Points technique combined with piecewise constancy zones (explained in section 2.2).
A value of 1×10−4 m−1 was assigned to the Specific Storage (Ss) over the entire model, while the Specific Yield (Sy) was assumed to be 20% for the permeable portions of the Aquifer Groups and 5% for the less permeable portions, which are values suitable for the lithology of the area.
2.3 Global sensitivity analysis, model calibration and validation
A sensitivity analysis was carried out to evaluate the influence of parameters and boundary conditions on model outputs. A global sensitivity analysis (GSA) was performed using the method of Morris [also called Elementary Effects Test (Almahawis et al., 2024; Abbas et al., 2025; Morris, 1991)], by means of the pestpp-sen tool, which is part of the PEST++ suite (White et al., 2020). The Morris method evaluates parameter sensitivity through two complementary indices. The first is the mean of the absolute elementary effects (μ*), following the modification proposed by Campolongo et al. (2005), which measures the overall influence of each input factor on the model output. The second is the standard deviation of the elementary effects (σ), which reflects the presence of nonlinear effects and/or interactions among input factors, with larger values indicating stronger nonlinearity or parameter interactions (Abbas et al., 2024).
The parameters considered in the GSA included CH heads (main rivers); GHB heads and conductance (northern and southern boundary conditions); drain stages and conductance (lowland springs); river stages and conductance (internal streams); irrigation canal losses; well pumping rates; recharge; specific yield; and hydraulic conductivity represented by Pilot Points (PPs), for a total of 618 adjustable parameters.
Hydraulic conductivity was parameterized using the Pilot Point (PP) method coupled with four piecewise-constancy zones to account for abrupt changes in sediment properties, such as clay lenses embedded within a more permeable matrix. Pilot Points are discrete locations within the model domain where hydraulic conductivity is estimated during calibration within bounds consistent with the geological setting. The estimated values are then interpolated by kriging separately within each piecewise-constancy zone, using only the Pilot Points belonging to the same zone. This approach produces a cell-by-cell hydraulic conductivity field while preserving sharp contrasts across hydrostratigraphic boundaries. The four piecewise-constancy zones were defined to represent the main hydrostratigraphic units and lithological contrasts. Zones 1 and 3 correspond to the permeable sediments of Aquifers A and B, respectively; Zone 2 represents the clay aquitard separating the two aquifers, whereas Zone 4 represents clay lenses embedded within Aquifer B. The same parameterization was adopted in Alberti et al. (2025) and is kept in the present study. A total of 607 PPs were distributed throughout the model domain and across the model layers to estimate the spatial distribution of hydraulic conductivity while reproducing the observed hydraulic heads.
Hydraulic conductivity (at the Pilot Points) and specific yield were varied directly using the absolute parameter values, whereas all remaining parameters were varied by applying multiplicative factors to their initial values. Input files were modified using a framework developed in FloPy (Bakker et al., 2016; Hughes et al., 2024), a Python package able to create, read and modify MODFLOW input and output files (the parameter list and their range of values are presented in Supplementary Table S2). To reduce the computational burden of the sensitivity analysis, the 618 adjustable parameters were grouped according to the model feature they represent. For example, all parameters defining the General Head Boundaries (GHB), including boundary heads and conductances, were assigned to a single group, while all river, drain, well, irrigation canal, recharge, storage, and hydraulic conductivity parameters were grouped in the same manner. Consequently, the Morris analysis evaluated the sensitivity of these nine model components rather than that of each individual parameter.
The model parameters were calibrated by means of a combination of Trial and Error (T&E) and assisted calibration by means of PEST_HP (Doherty, 2020). The parameters modified during the history-matching process were the recharge in urban areas (T&E) and the hydraulic conductivity (assisted calibration). The calibration employed Singular Value Decomposition (SVD) and Tikhonov regularization (Doherty et al., 2011; Fienen et al., 2009; Aissat et al., 2023).
The observations used in the calibration process consisted in a dataset of 3′510 hydraulic head values, measured in 297 monitoring points distributed within the model domain and among the different layers. The hydraulic head data was retrieved from the official database of ARPA Lombardia (Lombardy Environmental Protection Agency), Water Manager databases (MM S.p.a., CAP Holding S.p.a) and project-related field measurements.
After the calibration process, the model was validated on a different observation dataset, using head measurements available for the year 2019. The dataset was retrieved from ARPA Lombardia and MM S.p.a. groundwater level measurement databases, for a total of 712 observations measured in 134 monitoring wells. The model was then adapted to the 2019 conditions. The vertical recharge was estimated with the same methodology described above, and the pumping rates were updated.
2.4 Irrigation reduction scenarios
Once the transient model was implemented and calibrated, it was used to evaluate the groundwater system long term response to a decrease in the recharge component coming from irrigation. The model was then extended in time util 2050, cyclically replicating the conditions of the five-year period 2014–2018 in terms of meteoric recharge, well extractions and the other internal conditions. The repetition of the same meteoric recharge rates was chosen since the 2014–18 data vary from a wet period in 2014 to a relatively drier period in 2017, thus representing different recharge conditions over the years.
The scenarios consider a reduction of irrigation water mirroring possible changes in irrigation practices or a decrease in blue water availability, to assess their possible impact on groundwater levels; however, they are not meant to make actual predictions of groundwater levels but to assess the effects of such changes in a long-term perspective. Such changes were abruptly applied to the entire irrigated area after the first 5 years (i.e., after 2018), kept as in the calibrated model. In agricultural areas, three different reductions in irrigation recharge were evaluated: 15, 30, and 50% from initial values (scenarios Irr15, Irr30 and Irr50, respectively).
3 Results and discussion
3.1 Global sensitivity analysis
Morris’s sensitivity results on the considered parameters, presented in Figure 4, show that the most influencing parameter (i.e., the one with the largest value of ) is the hydraulic conductivity (HK), which also shows the highest value of , therefore suggesting nonlinearity effects or factor interactions linked to this parameter group. Other major influencing parameters are the storage (Specific Yield) and the internal stream parameters (stage and conductance), while the others show a lower importance. Recharge was varied using a multiplication factor within a ± 5% range (i.e., 0.95 to 1.05). The Morris sensitivity analysis indicated that recharge has a limited influence on the model outputs. Although recharge is generally considered an uncertain parameter, the availability of detailed information from the Irrigation Consortiums and distributed meteorological data used for its estimation reduced its uncertainty in the present study. Consequently, recharge was considered of secondary importance during calibration, which therefore focused primarily on the hydraulic conductivity parameters.
Figure 4
3.2 Transient model calibration and validation
The statistics obtained by comparing simulated and observed values are reported in Table 2, while the scatterplot between the observed and simulated data is shown in Figure 5.
Table 2
| Statistics | Calibration period (2014–2018) | Validation period (2019) |
|---|---|---|
| Residual mean [m] | 0.07 | 0.18 |
| Absolute residual mean [m] | 1.08 | 1.22 |
| Residual standard deviation [m] | 1.45 | 1.64 |
| RMSE [m] | 1.46 | 1.65 |
| Minimum residual [m] | −6.16 | −6.05 |
| Maximum residual [m] | 5.82 | 4.8 |
| Scaled absolute mean [%] | 1.0 | 1.0 |
| Scaled residual standard deviation [%] | 1.3 | 2.0 |
Statistics comparing observed and simulated data on the calibration and validation periods.
Figure 5
Considering the model extent (1880 km2) and its complexity, an absolute residual mean value of 1.08 m is deemed a satisfactory indicator of the model’s capability to reproduce groundwater levels and its variations over time. The model’ scaled residual standard deviation value (1.3%) is lower than the value commonly recognized as an indicator of good model reproduction (10%) (Anderson et al., 2015), thus indicating a good performance. The scatter plot shows a good match between the model simulation and the measurements used in the targets, without showing any trend departing from the bisector line. The mass balance error was lower than 1% in each SPs.
In Figure 6b, observed and simulated time series of four target points are shown (graphs for additional four points are shown in Supplementary Figure S2). An acceptable correspondence in terms of seasonal variations and long-term trend is observed in most of the analyzed series. Differences between the computed and observed were found. This could be due to several factors, such as uncertainties in monitoring data (most of these monitoring points are pumping wells usually active) and the presence of wells close to the monitoring points. The model is able to reproduce the observed behaviors although it slightly overestimated some wells while underestimating others inside Milan’s municipality boundaries.
Figure 6
Two examples of simulated piezometric heads for Aquifer Group A are shown in Figure 6a, in a period of high average piezometric head (December 2014) and low average piezometric head (December 2018). The computed heads relative to Aquifer Group B are shown in Supplementary Figure S3. The greatest differences in terms of heads are located within Milan, with differences reaching a maximum of 2 to 3 meters. At the edges of the domain no relevant differences are present, due to the steady-state nature of the boundary conditions. The calibrated hydraulic conductivity field is shown in Supplementary Figure S4.
The validation was carried out on observations measured in 2019 and showed a good agreement between observed and computed values, with an absolute residual mean of 1.22 m and a RMSE of 1.65 m (full statistics in Table 3 and calibration graph in Supplementary Figure S5), which are comparable to the calibration ones. Considering the model extension and complexity, the model can thus be considered validated.
Table 3
| Flow component | Recharge term | Average flux (m3/s) | Percent contribution |
|---|---|---|---|
| Inflow | Well* | 3.5 | 5% |
| River | 2.0 | 3% | |
| GHB | 14.3 | 23% | |
| CH** | 1.9 | 3% | |
| Storage | 10.9 | 17% | |
| Recharge | 30.9 | 49% | |
| Total | 63.4 | - | |
| Outflow | Well | 20.4 | 32% |
| River | 5.5 | 9% | |
| Drain | 0.7 | 1% | |
| GHB | 9.0 | 14% | |
| CH** | 18.3 | 29% | |
| Storage | 9.5 | 15% | |
| Total | 63.4 | - |
Model water budget in terms of average fluxes and percent contribution to the overall balance.
* This inflow represents the main irrigation canals leakages; ** Ticino and Adda rivers.
3.3 Mass balance analysis
A mass balance analysis was conducted on the calibrated model and on the recharge components to understand the relative importance of the different stresses applied to the model.
Figure 7 and Table 2 clearly show how the recharge is the most important inflow voice of the water budget (49% of the total water budget). The graph shows higher recharge rates starting from 2014 (1′646 mm of cumulated precipitation at the Milano Lambrate station, ID 2006), that progressively decrease due to less rainy years (e.g., only 685 mm of cumulated precipitation in 2017). However, it can be observed that the recharge increases every year from April to September due to irrigation. Overall, irrigation-derived recharge covered 57% of the total recharge volume provided, while 34% was provided by meteoric recharge and 8% by urban losses. To this volume, additional 3.5 m3/s is added by the main irrigation canals’ losses (WELL* in Table 2). The flood irrigation technique, the most used in the study area, plays thus a fundamental role in the aquifer recharge rates, as well as the canals losses themselves (Alberti et al., 2008). The relation between the irrigation recharge and the storage is clear. During the spring and summer time the storage component increases due to the water coming from the irrigation (negative fluxes since the direction is from the aquifer to the storage), while the aquifer releases water during the wintertime when the recharge is limited (positive fluxes, release from storage to the aquifer). The storage accounts for 17% of the total inflow, on average. The inflow from the upstream sector of the model, represented by GHBs, is almost constant during the simulation period and accounts for 23% of the total budget. Irrigation canal losses, the main rivers (Adda and Ticino, represented by CHs) and minor rivers are less important terms as inflow terms.
Figure 7
The most important outflow term for the area is groundwater abstraction by wells, which contributes 32% of the total. The main rivers (CH) and the southern boundary conditions (GHB) are the most important drainage elements in the domain, accounting for the 29 and 14% of the total outflow, respectively. The storage outflow is almost balanced with the inflow from storage (15%), while drains and minor rivers are less important in the overall budget.
The simulated flow rates of lowland springs are very small and vary during the year, reaching the maximum value in the irrigation season (3.5 m3/s) and the minimum one in winter, when many of them become dry. However, considering the preliminary results of the ongoing research of the MAURICE project (CE0100184) Colombo et al. (2024), the flow rate can be underestimated. This is mainly due to the lack of information regarding the lowland springs: spring elevations, spring-bed hydraulic conductivity are mainly retrieved from past studies (Bischetti et al., 2012) and flow rates data is scarce and present for only one point in time or for brief past time series. Nevertheless, the water budget should be only slightly influenced by the lowland springs discharge, since most of the water is used for field irrigation few hundred meters downstream. Furthermore, based on the sensitivity analysis results, lowland springs parameters do not significantly affect the computed heads.
The model reports a groundwater budget for aquifers A and B that aligns with previous findings for the Adda-Ticino basin (Alberti et al., 2016), thereby reinforcing the consistency of previous assessments and confirming that the recharge and withdrawal components are of the same order of magnitude as the incoming groundwater flow. The withdrawal-to-aquifer recharge ratio results equal to 0.66 for the period 2014–2018, suggesting a sustainable use of the underground water system.
3.4 Change in irrigation practices scenarios
Maps representing groundwater head difference were realized by subtracting 2018 heads from 2050 heads (head2050-head2018) and are shown in Figure 8. The first scenario is the most conservative one, considering a 15% reduction in recharge coming from irrigation from the year 2018 through 2050. Figure 8a shows that this results in a general lowering of groundwater levels over most of the model domain. The most affected areas are those located on the western and the eastern sides of the City of Milan, the zones with the largest agricultural areas. Head variations are present also in the North and South limits of the model, thanks to the flexibility to internal condition variations of the GHB conditions implemented at the boundaries. However, the lowering of the water table expands less towards the north boundary as there are no irrigated areas upstream of the Villoresi Canal (Figure 3b). The entire Milan municipality is subject to head variations between 1 m and 1.5 m. In the second scenario, an intermediate decrease in irrigation recharge, equal to 30%, was assumed. In Figure 8b it can be observed that the areas that undergo greater lowering are still those located on the sides of Milan, but in this case, the influence on groundwater levels within the city borders is greater. Compared to the previous scenario, it is evident how the presence of the Olona and Lambro rivers influences the hydraulic heads, moderating the lowering in the southernmost area of the domain. In both the northern and southern boundaries of the model, a maximum water table decrease of 1.8 m is recorded. In the city of Milan, the water table is subject to a decrease varying from 1.9 m to the East to almost 4 m to the western boundary.
Figure 8
The third scenario was implemented by assuming a drastic reduction in irrigation recharge, equal to 50% of that prior to 2018. The variation in terms of recharge can be observed in Supplementary Figure S6. Such a scenario, while unlikely in the immediate future, allows to identify the response of the system to an extreme situation. The differences between the final (2050) and initial (2018) model heads are shown in Figure 8c. Consistent changes in piezometric levels over most of the modeling domain are visible in this scenario. The feeding relationship that rivers have with the aquifer of interest is even more visible, and it’s clear the buffering effect of the Adda and Ticino rivers. The most affected areas (on the sides of the city of Milan) increase in extension and the decrease in groundwater levels reach values greater than 6 m, with peaks above 7 m. Inside the city of Milan, the lowering ranges from a minimum of 3 m (in the easternmost area) to a maximum of 6.7 m (in the westernmost portion).
Figure 9 shows groundwater level hydrographs at some selected locations within the model domain (Figure 9a). The hydrographs reflect what was already observed in the maps, showing some areas which are more influenced by the recharge rate reduction. As expected, agricultural areas are the most sensitive to irrigation practices change. For example, comparing the Base and IRR50 scenarios, in monitoring point o106 shows a groundwater level decline from average values of 133 m a.s.l. to 127 m a.s.l., with the largest drawdown occurring in the first 10 years. In contrast, for the same scenarios, a smaller response is simulated in urbanized areas. At monitoring point o103, the average groundwater level decreases from about 133 m a.s.l. to approximately 131.5 m a.s.l., with a more gradual decline between 2018 and 2032. In general, the graphs show a stabilization of the levels in each scenario, which is reached earlier by the less critical scenarios (e.g., Irr15). Therefore, results suggest that the groundwater table will reach stability after some years, but with lower groundwater levels with a potential drawdown higher than 7 m at some locations in the Irr50 scenario. This could likely affect shallow wells, dry-up wetlands and most of the lowland springs and increase the drainage from surface water bodies.
Figure 9
4 Conclusion
The development of numerical models to simulate groundwater flow at the regional scale provides tools to support water managers in the quantitative management of groundwater resources. This transient state model is the first of its kind at the Adda-Ticino basin scale and it is a fundamental tool to quantitatively assess the effects on groundwater by the variation of climatic and anthropogenic forcings within the territory, such as the effects of climate change (increased precipitation or drought periods), the change in land use, the implementation of new infrastructure for the extraction of water for drinking purposes or changes in irrigation practices. In this study, the impacts of potential future optimization procedures of irrigation systems may have on groundwater levels below the city of Milan were evaluated.
The calibration results show that the model is able to reproduce with good accuracy the evolution of groundwater levels over time caused by forcing variables (recharge and withdrawal). Furthermore, the validation procedure carried on a different observation dataset showed satisfactory reliability in reproducing other periods. Therefore, the model can be considered suitable to perform simulation of future scenarios and, furthermore, to be used as a tool to support the decisions processes of water managers at basin scale and the development of water safety plans.
The water budget analysis showed the sustainability of groundwater use in the area with an increase of the downstream water flows compared to the incoming flows. The implemented scenarios of changes in irrigation practices showed their potential effects, showing a significant and dangerous reduction in the overall groundwater levels. The study highlights the importance of traditional irrigation methodologies in sustaining the groundwater level in the Milan area. Although surface water saving policies should be evaluated, before their adoption it is necessary to consider the potential effects on groundwater availability as well, in this case a loss of a significant fraction of its recharge. Furthermore, adopting sprinkler or drip irrigation systems could also lead to the loss of the capillary irrigation canals’ inherent heritage and an increase in the number of wells used for irrigation, escalating the stresses on the underground water resources.
This work represents a starting point to assess future piezometric levels in the Milan area considering hydrological balance changes (e.g., climate change effects), to develop robust mitigation and adaptation measures to ensure water availability at the basin scale.
Future research and improvements in the hydrogeological characterization of the study area will further enhance the reliability of the proposed model and reduce the uncertainty of its predictions. Once data covering the dry period and the subsequent wet period (2022–2025) become available, the model can be validated under these extreme hydrological conditions, thereby increasing confidence in its predictive capability. In addition, pumping test data held by Public Authorities, which are currently not readily available, would help constrain the hydraulic conductivity distribution and reduce parameter uncertainty. Additional information on the hydraulic characteristics of the lowland springs would improve the representation of groundwater discharge processes, while a more detailed characterization of the effective porosity distribution would improve the simulation of groundwater responses to transient hydrological forcing. Finally, a more accurate description of irrigation water distribution and irrigation practices would allow a better quantification of irrigation recharge, further reducing its uncertainty.
Statements
Data availability statement
The data should be asked to the competent authority, wells positions and abstractions to Regione Lombardia (www.regione.lombardia.it), irrigation data to the Consorzio Est Ticino Villoresi (www.etvilloresi.it) and Consorzio Bonifica Muzza Bassa Lodigiana (www.muzza.it), while piezometric and hydrometeorological data can be retrieved from the ARPA Lombardia website (www.arpalombardia.it) or asked to the water managers (MM (www.mmspa.eu) and CAP (www.gruppocap.it)).
Author contributions
PM: Software, Methodology, Conceptualization, Writing – review & editing, Validation, Formal analysis, Writing – original draft. PC: Writing – review & editing, Data curation, Methodology, Writing – original draft, Software, Visualization. LC: Writing – review & editing, Methodology. LA: Funding acquisition, Validation, Supervision, Writing – review & editing, Conceptualization.
Funding
The author(s) declared that financial support was not received for this work and/or its publication.
Acknowledgments
The authors thank MM S.p.a., CAP Holding S.p.a., ARPA Lombardia, Consorzio di Bonifica Est Ticino Villoresi and Regione Lombardia for their collaboration in providing data.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that Generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frwa.2026.1872897/full#supplementary-material
References
1
AbbasS.BaileyR.WhiteJ. T.ArnoldJ.WhiteM. (2025). Estimation of groundwater storage loss using surface–subsurface hydrologic modeling in an irrigated agricultural region. Sci. Rep.15, 1–24. doi: 10.1038/s41598-025-92987-6
2
AbbasS. A.BaileyR. T.WhiteJ. T.ArnoldJ. G.WhiteM. J.ČerkasovaN.et al. (2024). A framework for parameter estimation, sensitivity analysis, and uncertainty analysis for holistic hydrologic modeling using SWAT+. Hydrol. Earth Syst. Sci.28, 21–48. doi: 10.5194/hess-28-21-2024
3
AbbottB. W.BishopK.ZarnetskeJ. P.HannahD. M.FreiR. J.MinaudoC.et al. (2019). A water cycle for the Anthropocene. Hydrol. Process.33, 3046–3052. doi: 10.1002/hyp.13544
4
AissatR.PryetA.SaltelM.DupuyA. (2023). Comparison of different pilot point parameterization strategies when measurements are unevenly distributed in space. Hydrogeol. J.31, 2381–2400. doi: 10.1007/s10040-023-02737-z
5
AlbertiL.AntelmiM.ObertoG.La LicataI.MazzonP. (2022). Evaluation of fresh groundwater Lens volume and its possible use in Nauru Island. Water14:20. doi: 10.3390/w14203201,
6
AlbertiL.CantoneM.ColomboL.LombiS.PianaA. (2016). Numerical modeling of regional groundwater flow in the adda-Ticino Basin: advances and new results. Rend. Online Soc. Geol. Ital.41/2016, 10–13. doi: 10.3301/ROL.2016.80
7
AlbertiL.ColomboL.FormentinG. (2018). Null-space Monte Carlo particle tracking to assess groundwater PCE (tetrachloroethene) diffuse pollution in north-eastern Milan functional urban area. Sci. Total Environ.621, 326–339. doi: 10.1016/j.scitotenv.2017.11.253,
8
AlbertiL.FrancaniV. (2001). Studio idrogeologico sulle cause del sollevamento della falda nell’area milanese. Geoing. Ambient. Miner.104, 257–264. Available online at: https://re.public.polimi.it/handle/11311/247013
9
AlbertiL.FrancaniV.InnocentiI. (2008). La Gestione Integrata Delle Risorse Idriche Lombarde: Valutazione Del Contributo Alla Ricarica Da Parte Della Rete Irrigua Una Nuova Geologia per La Lombardia Milano
10
AlbertiL.FrancaniV.MasettiM.ParriA. (2000). Valutazione del livello massimo raggiungibile dalla falda nel Comune di Milano. Quad. Geol. Appl.2, 13–28. Available online at: https://re.public.polimi.it/handle/11311/247013
11
AlbertiL.MazzonP.ColomboL.CantoneM.AntelmiM.MarelliF.et al. (2025). Enhancing groundwater resource Management in the Milan Urban Area through a robust stratigraphic framework and numerical Modeling. Water17:2. doi: 10.3390/w17020165,
12
AliA.HussainT.ZahidA. (2025). Smart irrigation technologies and prospects for enhancing water use efficiency for sustainable agriculture. AgriEngineering7:106. doi: 10.3390/agriengineering7040106
13
AlmahawisM. K.BaileyR. T.AbbasS. A.ArnoldJ. G.WhiteM. J. (2024). Investigating the impact of irrigation practices on hydrologic fluxes in a highly managed River Basin. Agric. Water Manag.301:108954. doi: 10.1016/j.agwat.2024.108954,
14
AndersonM. P.WoessnerW. W.HuntR. J. (2015). Applied Groundwater Modeling: Simulation of Flow and Advective Transport. Academic press. Available online at: https://books.google.com/books?hl=it&lr=&id=jcOcBAAAQBAJ&oi=fnd&pg=PP1&dq=anderson+and+woessner&ots=a2NYKcNpaK&sig=f-1LLN0rhvpqFskdoLmuPxi1mwM (Accessed August 05, 2024)
15
AntelmiM.AlbertiL.BarbieriS.PandayS. (2021a). Simulation of thermal perturbation in groundwater caused by borehole heat exchangers using an adapted CLN package of MODFLOW-USG. J. Hydrol.596:126106. doi: 10.1016/j.jhydrol.2021.126106,
16
AntelmiM.MazzonP.HöhenerP.MarchesiM.AlbertiL. (2021b). Evaluation of MNA in a chlorinated solvents-contaminated aquifer using reactive transport Modeling coupled with isotopic fractionation analysis. Water13:21. doi: 10.3390/w13212945,
17
AntelmiM.RenoldiF.AlbertiL. (2020). Analytical and numerical methods for a preliminary assessment of the remediation time of pump and treat systems. Water12:10. doi: 10.3390/w12102850,
18
BabuR.ParkN.YoonS.KulaT. (2018). Sharp Interface approach for regional and well scale Modeling of Small Island freshwater Lens: Tongatapu Island. Water10:11. doi: 10.3390/w10111636,
19
BakerE. A.CappatoA.TodeschiniS.TamelliniL.SangalliG.RealiA.et al. (2022). Combining the Morris method and multiple error metrics to assess aquifer characteristics and recharge in the lower Ticino Basin, in Italy. J. Hydrol.614:128536. doi: 10.1016/j.jhydrol.2022.128536,
20
BakkerM.PostV.LangevinC. D.HughesJ. D.WhiteJ. T.StarnJ. J.et al. (2016). Scripting MODFLOW model development using Python and FloPy. Groundwater54, 733–739. doi: 10.1111/gwat.12413,
21
BarbieriS.AntelmiM.MazzonP.RizzoS.AlbertiL. (2025). Numerical Modeling of the groundwater temperature variation generated by a ground-source heat pump system in Milan. Appl. Sci.15:10. doi: 10.3390/app15105522,
22
BecciuG.PaolettiA.. (2010). Fondamenti di Costruzioni Idrauliche. ITA. Available online at: https://re.public.polimi.it/handle/11311/562128
23
BerettaG. P.AvanziniM.PagottoA. (2004). Managing groundwater rise: experimental results and modelling of water pumping from a quarry lake in Milan urban area (Italy). Environ. Geol.45, 600–608. doi: 10.1007/s00254-003-0918-7
24
BischettiG. B.FumagalliN.PiantanidaE. V.et al. (2012). Tutela e valorizzazione dei fontanili del territorio lombardo—FonTe, vol. 144 Quaderni della Ricerca. Milan, Italy: Regione Lombardia.
25
BonomiT.Del RossoF.FumagalliL.CanepaP. (2010). Assessment of groundwater availability in the Milan province aquifers. Mem. Descr. Carta Geol. d’It, 31–40. Available online at: https://www.isprambiente.gov.it/files/pubblicazioni/periodicitecnici/memorie/memoriexc/90-bonomi.pdf
26
BossardM.FeranecJ.OťaheľJ.. (2000). Corine Land Cover—Technical Guide. Copenhagen, Denmark: European Environment Agency.
27
BrakefieldL. K.WhiteJ. T.HoustonN. A.ThomasJ. V. (2015). “Updated numerical model with uncertainty assessment of 1950-56 drought conditions on brackish-water movement within the Edwards aquifer, San Antonio, Texas,” in Scientific Investigations Report, Nos, (San Antonio, Texas: U.S. Geological Survey).
28
BrouwerC.PrinsK.KayM.HeibloemM. (1989). Irrigation Water Management: Irrigation Scheduling. FAO Irrigation and Drainage24. Available online at: https://www.fao.org/4/t7202e/t7202e00.htm
29
BultéM.DurenT.BouhonO.PetitclercE.AgnielM.DassarguesA. (2021). Numerical Modeling of the interference of thermally unbalanced aquifer thermal energy storage Systems in Brussels (Belgium). Energies14:19. doi: 10.3390/en14196241,
30
CampolongoF.CariboniJ.SchoutensW.. (2005). “Enhancing the Morris method.” Available online at: https://publications.jrc.ec.europa.eu/repository/handle/JRC31319 (Accessed June 23, 2025)
31
ChenC.-S.TuC.-H.ChenS.-J.ChenC.-C. (2016). Simulation of groundwater contaminant transport at a decommissioned landfill site—a case study, Tainan City, Taiwan. Int. J. Environ. Res. Public Health13. doi: 10.3390/ijerph13050467,
32
Città Metropolitana di Milano (Dataset della categoria AMBIENTE). (2024) “Isopieze I falda 2007-2013 (Isopiezometric maps from 2007 to 2013).” Available online at: https://www.cittametropolitana.mi.it/open_data/dataset-pubblicati-per-categoria/ambiente/index.html (Accessed August 26, 2024)
33
ColomboL.AlbertiL.MazzonP.AntelmiM. (2020). Null-space Monte Carlo particle backtracking to identify groundwater tetrachloroethylene sources. Front. Environ. Sci.8. doi: 10.3389/fenvs.2020.00142
34
ColomboL.AlbertiL.MazzonP.FormentinG. (2019). Transient flow and transport modelling of an historical CHC source in north-West Milano. Water11:9. doi: 10.3390/w11091745,
35
ColomboL.GattinoniP.ScesiL. (2018). Stochastic modelling of groundwater flow for Hazard assessment along the underground infrastructures in Milan (northern Italy). Tunn. Undergr. Space Technol.79, 110–120. doi: 10.1016/J.TUST.2018.05.007,
36
ColomboP.MazzonP.AlbertiL. (2024). Off-season irrigation as a climate adaptation strategy for future groundwater management in northern Italy. Adv. Geosci.64, 27–31. doi: 10.5194/adgeo-64-27-2024
37
CornesR. C.van der SchrierG.van den BesselaarE. J. M.JonesP. D. (2018). An ensemble version of the E-OBS temperature and precipitation data sets. J. Geophys. Res. Atmos.123, 9391–9409. doi: 10.1029/2017JD028200
38
CoulonC.LemieuxJ.-M.PryetA.BayerP.YoungN. L.MolsonJ. (2022). Pumping optimization under uncertainty in an island freshwater lens using a sharp-interface seawater intrusion model. Water Resour. Res.58:e2021WR031793. doi: 10.1029/2021WR031793
39
CoulonC.WhiteJ. T.PryetA.GatelL.LemieuxJ.-M. (2024). An ensemble-based approach for pumping optimization in an island aquifer considering parameter, observation and climate uncertainty. Hydrol. Earth Syst. Sci.28, 303–319. doi: 10.5194/hess-28-303-2024
40
DahlkeH. E.LaHueG. T.MautnerM. R. L.et al. (2018). “Chapter eight—managed aquifer recharge as a tool to enhance sustainable groundwater Management in California: examples from field and Modeling studies,” in Advances in Chemical Pollution, Environmental Management and Protection, eds. FriesenJ.Rodríguez-SinobasL., vol. 3 (Advanced Tools for Integrated Water Resources Management. Amsterdam, Netherland: Elsevier).
41
De CaroM.PericoR.CrostaG. B.FrattiniP.VolpiG. (2020). A regional-scale conceptual and numerical groundwater flow model in fluvio-glacial sediments for the Milan metropolitan area (northern Italy). J. Hydrol. Reg. Stud.29:100683. doi: 10.1016/j.ejrh.2020.100683
42
De PascaleS.CostaL. D.ValloneS.BarbieriG.MaggioA. (2011). Increasing Water Use Efficiency in Vegetable Crop Production: From Plant to Irrigation Systems Efficiency. HortTechnology21:301–308.
43
DelfiniI.ChahoudA.ZamrskyD.MontanariA. (2024). “Assessment of Climate Change Impacts on Groundwater Exploitation Sustainability in the Emilia-Romagna Region (Italy).” SSRN Scholarly Paper No. 4918779. Social Science Research Network. Available online at: https://papers.ssrn.com/abstract=4918779 (Accessed January 01, 2025)
44
DelsmanJ. R.MulderT.Romero VerasteguiB.BootsmaH.ZitmanP.HuizerS.et al. (2023). Reproducible construction of a high-resolution national variable-density groundwater salinity model for the Netherlands. Environ. Model. Softw.164:105683. doi: 10.1016/j.envsoft.2023.105683
45
DohertyJ.. (2020). PEST—Model-Independent Parameter Estimation Manual. Watermark Numerical Computing. Available online at: https://www.epa.gov/hydrowq/pest-model-independent-parameter-estimation-manual-4th-ed (Accessed August 26, 2024)
46
DohertyJ. E.FienenM. N.HuntR. J. (2011). “Approaches to highly parameterized inversion: pilot-point theory, guidelines, and research directions,” in Scientific Investigations Report, (USA: U.S. Geological Survey).
47
DyringM.HofmannH.StantonD.MossP.FroendR. (2023). Ecohydrology of coastal aquifers in humid environments and implications of a drying climate. Ecohydrology16:e2491. doi: 10.1002/eco.2491
48
European Commission. (2026). Water—Agriculture and Rural Development—European Commission. Available online at: https://agriculture.ec.europa.eu/cap-my-country/sustainability/environmental-sustainability/natural-resources/water_en
49
European Union. Regulation (EU) 2020/741, 177 OJ L (2020). Available online at: http://data.europa.eu/eli/reg/2020/741/oj (Accessed June 18, 2026)
50
FeinsteinD. T.HuntR. J.ReevesH. W. (2010). “Regional groundwater-flow model of the Lake Michigan Basin in support of Great Lakes Basin water availability and use studies,” in Scientific Investigations Report 2010–5109, (USA: U.S. Geological Survey).
51
FienenM. N.MuffelsC. T.HuntR. J. (2009). On constraining pilot point calibration with regularization in PEST. Groundwater47, 835–844. doi: 10.1111/j.1745-6584.2009.00579.x,
52
FreezeR. A.WitherspoonP. A. (1966). Theoretical analysis of regional groundwater flow: 1. Analytical and numerical solutions to the mathematical model. Water Resour. Res.2, 641–656. doi: 10.1029/WR002i004p00641
53
GattinoniP.ScesiL. (2017). The groundwater rise in the urban area of Milan (Italy) and its interactions with underground structures and infrastructures. Tunn. Undergr. Space Technol.62, 103–114. doi: 10.1016/j.tust.2016.12.001,
54
GiudiciM.FogliaL.ParraviciniG.PonziniG.SincichB. (2000). A quasi three dimensional model of water flow in the subsurface of Milano (Italy): the stationary flow. Hydrol. Earth Syst. Sci.4, 113–124. doi: 10.5194/hess-4-113-2000
55
GiulianoG. (1995). Ground water in the Po Basin: some problems relating to its use and protection. Sci. Total Environ. Toxicol.171, 17–27. doi: 10.1016/0048-9697(95)04682-1
56
GorlaM. (2001). Acquiferi Confinati Profondi della Provincia di Milano: Sviluppi delle Conoscenze Idrogeologiche e Idrochimiche ai fini di un loro uso Sostenibile. Deep Confined Aquifers within the Province of Milan: Hydrogeological and Hydrochemical Evaluation for their Sustainable use. Part 2nd–Geochemical Features, Pollutants and Conclusions.
57
GorlaM.SimonettiR.RighettiC. (2016). Basin-scale hydrogeological, geophysical, geochemical and isotopic characterization: an essential tool for building a decision support system for the sustainable management of alluvial aquifer systems within the provinces of Milan and Monza-Brianza. Acque Sotterranee Ital. J. Groundwater5:3. doi: 10.7343/as-2016-213
58
HarbaughA. W. (2005). MODFLOW-2005, the US Geological Survey Modular Ground-Water Model: The Ground-Water Flow Process Vol. 6 US Department of the Interior, US Geological Survey Reston, VA, USA. Available online at: http://wwwbrr.cr.usgs.gov/hill_tiedeman_book/documentation/MODFLOW-MODPATH-ModelViewer/MF2005-tma6a16.pdf
59
HargreavesG. H.SamaniZ. A. (1985). Reference crop evapotranspiration from temperature. Appl. Eng. Agric.1, 96–99. doi: 10.13031/2013.26773
60
HowellT. (2003). Irrigation Efficiency Encyclopedia of Water Science. Available online at: https://www.researchgate.net/publication/43256707_Irrigation_Efficiency (Accessed August 26, 2024)
61
HughesJ. D.LangevinC. D.PaulinskiS. R.LarsenJ. D.BrakenhoffD. (2024). “FloPy Workflows for Creating Structured and Unstructured MODFLOW Models,” in Groundwater, Python, v. 3.10.0.dev3. doi: 10.1111/gwat.13327
62
IslamM. B.FirozA. B. M.FogliaL.MarandiA.KhanA. R.SchüthC.et al. (2017). A regional groundwater-flow model for sustainable groundwater-resource management in the south Asian megacity of Dhaka, Bangladesh. Hydrogeol. J.25, 617–637. doi: 10.1007/s10040-016-1526-4
63
Istat (2014). 6° Censimento Generale Dell’Agricoltura—Utilizzo Della Risorsa Idrica a Fini Irrigui in Agricoltura. Available online at: https://www.istat.it/it/files/2014/11/Utilizzo_risorsa_idrica.pdf (Accessed July 16, 2024)
64
Istat. (2023). Profili delle Città Metropolitane. Available online at: https://www.istat.it/it/files/2023/02/Statistica-Focus-Citt%C3%A0-Metropolitane.pdf
65
IzukaS. K.RotzollK.NishikawaT. (2021). Volcanic Aquifers of Hawai‘i—Construction and Calibration of Numerical Models for Assessing Groundwater Availability on Kaua‘i, O‘ahu, and Maui Scientific Investigations Report, Nos. 2020–5126. Reston, VA: U.S. Geological Survey.
66
JacksonC. R.MeisterR.PrudhommeC. (2011). Modelling the effects of climate change and its uncertainty on UK chalk groundwater resources from an Ensemble of Global Climate Model Projections. J. Hydrol.399, 12–28. doi: 10.1016/j.jhydrol.2010.12.028,
67
JasechkoS.SeyboldH.PerroneD.FanY.ShamsudduhaM.TaylorR. G.et al. (2024). Rapid groundwater decline and some cases of recovery in aquifers globally. Nature625, 715–721. doi: 10.1038/s41586-023-06879-8,
68
LasagnaM.ManciniS.De LucaD. A. (2020). Groundwater hydrodynamic behaviours based on water table levels to identify natural and anthropic controlling factors in the Piedmont plain (Italy). Sci. Total Environ.716:137051. doi: 10.1016/j.scitotenv.2020.137051
69
LottiF.BorsiI.GuastaldiE.BarbagliA.BasileP.FavaroL.et al. (2021). Numerically enhanced conceptual modelling (NECoM) applied to the Malta mean sea level aquifer. Hydrogeol. J.29, 1517–1537. doi: 10.1007/s10040-021-02330-2
70
MasseroniD.GangiF.GalliA.CerianiR.De GaetaniC.GandolfiC. (2022). Behind the efficiency of border irrigation: lesson learned in northern Italy. Agric. Water Manag.269:107717. doi: 10.1016/j.agwat.2022.107717,
71
MorrisM. D. (1991). Factorial sampling plans for preliminary computational experiments. Technometrics33, 161–174. doi: 10.1080/00401706.1991.10484804
72
NiswongerR. G.PandayS.IbarakiM. (2011). MODFLOW-NWT, a Newton formulation for MODFLOW-2005. U. S. Geol. Surv. Tech. Methods6:44. Available online at: https://pubs.usgs.gov/tm/tm6a37/
73
OECD (2012). Redefining “Urban”: A New Way to Measure Metropolitan Areas. France, Paris: OECD.
74
OmagbonJ.DohertyJ.YehA.ColinaR.O'SullivanJ.McDowellJ.et al. (2021). Case studies of predictive uncertainty quantification for geothermal models. Geothermics97:102263. doi: 10.1016/j.geothermics.2021.102263
75
PerdikakiM.MakropoulosC.KalliorasA. (2022). Participatory groundwater modeling for managed aquifer recharge as a tool for water resources management of a coastal aquifer in Greece. Hydrogeol. J.30, 37–58. doi: 10.1007/s10040-021-02427-8
76
PiggotA. R. (2002). Regional groundwater Modeling within the Grand River watershed. NWRI contribution; 98-147. Environent Canada water Science and Technology Directorate. Available online at: https://publications.gc.ca/site/eng/9.872514/publication.html
77
PollicinoL. C.ColomboL.FormentinG.AlbertiL. (2021a). Stochastic modelling of solute mass discharge to identify potential source zones of groundwater diffuse pollution. Water Res.200:117240. doi: 10.1016/j.watres.2021.117240,
78
PollicinoL. C.MasettiM.StevenazziS.CristaldiA.RighettiC.GorlaM. (2021b). Multi-aquifer susceptibility analyses for supporting groundwater Management in Urban Areas. J. Contam. Hydrol.238:103774. doi: 10.1016/j.jconhyd.2021.103774,
79
PreviatiA.CrostaG. B. (2021). Characterization of the subsurface urban heat island and its sources in the Milan city area, Italy. Hydrogeol. J.29, 2487–2500. doi: 10.1007/s10040-021-02387-z
80
RedaelliA.BonomiT.SartiranaD.SinatraG.FeinsteinD. T.HuntR. J.et al. (2026). Changes in irrigation practices may deplete aquifers faster and more severely than meteorological droughts: a numerical Modeling approach. J. Hydrol.672:135337. doi: 10.1016/j.jhydrol.2026.135337,
81
Regione Lombardia. (2018). Dusaf 6.0—Uso Del Suolo 2018 | Open Data Regione Lombardia. Available online at: https://www.dati.lombardia.it/Territorio/Dusaf-6-0-Uso-del-suolo-2018/7rae-fng6/about_data (Accessed May 08, 2024)
82
Regione Lombardia (Geoportale della Lombardia). (2014). Piezometrie 2014 degli acquiferi superficiali e profondi. Available online at: https://www.geoportale.regione.lombardia.it/metadati (Accessed August 26, 2024)
83
Regione Lombardia and ENI-AGIP (2002). Geologia Degli Acquiferi Padani Della Regione Lombardia. Florence, Italy: S.EL.CA.
84
RinglebJ.SallweyJ.StefanC. (2016). Assessment of managed aquifer recharge through Modeling—a review. Water8:12. doi: 10.3390/w8120579,
85
RotirotiM.BonomiT.SacchiE.McArthurJ. M.StefaniaG. A.ZanottiC.et al. (2019). The effects of irrigation on groundwater quality and quantity in a human-modified hydro-system: the Oglio River basin, Po plain, northern Italy. Sci. Total Environ.672, 342–356. doi: 10.1016/j.scitotenv.2019.03.427,
86
RotirotiM.SacchiE.CaschettoM.ZanottiC.FumagalliL.BiasibettiM.et al. (2023). Groundwater and surface water nitrate pollution in an intensively irrigated system: sources, dynamics and adaptation to climate change. J. Hydrol.623:129868. doi: 10.1016/j.jhydrol.2023.129868,
87
RussoT. A.FisherA. T.LockwoodB. S. (2015). Assessment of managed aquifer recharge site suitability using a GIS and Modeling. Groundwater53, 389–400. doi: 10.1111/gwat.12213,
88
SarmaR.SinghS. K. (2022). A comparative study of data-driven models for groundwater level forecasting. Water Resour. Manag.36, 2741–2756. doi: 10.1007/S11269-022-03173-6/FIGURES/5
89
SartiranaD.ZanottiC.RotirotiM.CaschettoM.RedaelliA.BrunoS.et al. (2024). Urban water Management in Milan Metropolitan Area, a review. Acque Sotterranee Ital. J. Groundw.13:3. doi: 10.7343/as-2024-763,
90
SartiranaD.ZanottiC.RotirotiM.de AmicisM.CaschettoM.RedaelliA.et al. (2022). Quantifying groundwater infiltrations into Subway lines and underground Car parks using MODFLOW-USG. Water14:24. doi: 10.3390/w14244130,
91
TrefryM. G.MuffelsC. (2007). FEFLOW: a finite-element ground water flow and transport modeling tool. Groundwater45, 525–528. doi: 10.1111/j.1745-6584.2007.00358.x
92
United Nations (2022). The United Nations World Water Development Report 2022: Groundwater: Making the Invisible Visible. Paris, France: UNESCO.
93
VassenaC.RienznerM.PonziniG.GiudiciM.GandolfiC.DuranteC.et al. (2012). Modeling water resources of a highly irrigated alluvial plain (Italy): calibrating soil and groundwater models. Hydrogeol. J.20, 449–467. doi: 10.1007/s10040-011-0822-2
94
WestenbroekS. M.EngottJ. A.KelsonV. A.HuntR. J. (2018). SWB Version 2.0-a Soil-Water-Balance code for Estimating Net Infiltration and Other Water-Budget Components. Reston, VA:U.S. Geological Survey.
95
WhiteJ. T.HuntR. J.FienenM. N.DohertyJ. E. (2020). “Approaches to Highly Parameterized Inversion: PEST++ Version 5, a Software Suite for Parameter Estimation, Uncertainty Analysis, Management Optimization and Sensitivity Analysis,” in Techniques and Methods, vol. 7-C26 (Reston, VA: U.S. Geological Survey).
96
WoodwardD. E.HawkinsR. H.JiangR.HjelmfeltJ. A. V. M.QuanD. (2003). Runoff Curve Number Method: Examination of the Initial Abstraction Ratio, 1–10. doi: 10.1061/40685(2003)308
97
XuZ.SerataR.WainwrightH.DenhamM.MolinsS.Gonzalez-RaymatH.et al. (2022). Reactive transport modeling for supporting climate resilience at groundwater contamination sites. Hydrol. Earth Syst. Sci.26, 755–773. doi: 10.5194/hess-26-755-2022
98
YangL.ZhengC.AndrewsC. B.WangC. (2021). Applying a regional transport modeling framework to manage nitrate contamination of groundwater. Groundwater59, 292–307. doi: 10.1111/gwat.13047,
99
ZhangL.LiangL.ZhuJ.GanT.YangY.WangD. (2024). Impact of hydrogeological characteristics on the development of underground space resources in valley cities. Geofluids1:1052700. doi: 10.1155/2024/1052700
100
ZhouY.LiW. (2011). A review of regional groundwater flow modeling. Geosci. Front.2, 205–214. doi: 10.1016/j.gsf.2011.03.003
101
ZucaroR.CorapiA. (2009). Rapporto Sullo Stato Dell’irrigazione in Lombardia. Monitoraggio Dei Sistemi Irrigui Delle Regioni Centro Settentrionali. INEA. Available online at: https://sigrian.crea.gov.it/wp-content/uploads/2020/11/2009-Rapporto-sullo-stato-dellirrigazione-in-Lombardia.pdf
Summary
Keywords
groundwater management, groundwater modeling, irrigation management, MODFLOW, Northern Italy
Citation
Mazzon P, Colombo P, Colombo L and Alberti L (2026) Quantifying irrigation-driven recharge in a complex urban-agricultural system: the case of Milan, Italy. Front. Water 8:1872897. doi: 10.3389/frwa.2026.1872897
Received
05 May 2026
Revised
08 July 2026
Accepted
20 July 2026
Published
10 September 2026
Volume
8 - 2026
Edited by
Pragnaditya Malakar, Indian Institute of Science Education and Research Berhampur (IISER), India
Reviewed by
Aatish Anshuman, Indian Institute of Technology Bombay, India
Ilaria Delfini, University of Bologna, Italy
Updates
Copyright
© 2026 Mazzon, Colombo, Colombo and Alberti.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Paolo Colombo, paolo1.colombo@polimi.it
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.