Numerical modelling of land subsidence related to groundwater withdrawal in the Firenze-Prato-Pistoia basin (central Italy)

The Tuscany region of Italy is widely affected by subsidence, landslides and floods, which severely impact buildings and infrastructure. In particular, Firenze-Prato-Pistoia basin has a long experience of ground deformation related to groundwater withdrawal. European remote-sensing satellite (ERS) data collected since 1992 have revealed the presence of several subsiding areas in the basin such as the south-eastern portion of the city of Pistoia. Sentinel-1 persistent scatterer interferometry (PSI) measurements for 2015–2018 confirmed the long-term subsidence of this area, associated with intense horticulture (plant nurseries). At the same time, Sentinel-1 data revealed the unexpected movement of Pistoia historic center, which has always been considered stable in the past. To identify the complex relationship between aquifer conditions and ground displacement, a hydrogeologic model of the Pistoia aquifers was developed, applying an integrated modelling procedure. Hydrodynamic-parameter distributions, calibrated and validated by means of Sentinel-1 PSI measurements, suggest that subsidence in Pistoia area is probably related with the combined impacts of groundwater extraction and highly compressible aquitards. To evaluate the potential evolution of ground displacement, numerical simulations were extended until 2050, using regional and global climate model data, analyzing three different pumping-rate scenarios. This led to the development of several subsidence hazard maps of the city of Pistoia that display the influence of groundwater extraction in controlling land subsidence in the area. This study emphasizes the importance of developing proper groundwater management policies, especially in alluvial aquifers made of fine compressible sediments, in order to sustainably utilize underground freshwater resources and to avoid related side effects.


Introduction
Groundwater is a fundamental resource for freshwater supplies in many urban and rural areas all over the world, as it guarantees support for the agricultural economy and the availability of large amounts of drinkable water. Historically, groundwater was extracted only considering the human necessities, often with little regard for the effects of over pumping on natural aquifer systems. High groundwater extraction rates, greater than average aquifer recharge rates, are not sustainable for the environment in the long term, leading to the occurrence of many economic and environmental consequences such as water quality deterioration, increasing costs for replacement of groundwater pumps and/or wells, drainage of nearby rivers, and land subsidence (Fienen and Arshad 2016).
When groundwater pumping exceeds aquifer recharge for long periods, a decline of water levels usually occurs, with a relative decrease in pore-water pressure and an increase in effective stress. This phenomenon leads to the compression of aquifer-system materials and to gradual land-surface subsidence (Galloway and Burbey 2011). In San Joaquin valley, California (USA), uncontrolled groundwater extraction in the first decades of the last century, and poor management since then, has led to land subsidence by several meters, with a maximum measured ground displacement at the beginning of the 1970s of more than 8.5 m (Poland et al. 1975). In Mexico City, during the 20th century, the compaction of lacustrine deposits due to groundwater withdrawal caused several meters of land-surface deformation, with direct repercussions on human activities and aquifer vulnerability (Figueroa-Vega 1984;Hernández-Espriú et al. 2014).
Some of the most expansive consequences of land subsidence include damage to buildings and infrastructure, as well as the extension of coastal and inland flooding areas (Di Paola et al. 2017;Alberico et al. 2012). Buildings and underground facilities are under higher risks of damage if they are placed in locations characterized by important differential ground deformation (Shi et al. 2019).
To prevent major problems, the development of adequate groundwater management strategies should be one of the main priorities for local authorities and policy makers. The use of n u m e r i c a l g r o u n d w a t e r m o d e l s t o i d e n t i f y t h e hydrogeological condition of an area, together with interferometric synthetic aperture radar (InSAR) data that detect ground-surface displacement, provides a powerful tool to analyze aquifer threats and these methods have been extensively covered in literature in recent years (Ezquerro et al. 2018;Arjomandi et al. 2018). However, one of the biggest uncertainties found in the proposed methodologies is represented by the traditional approach of water resource modelling. Such modelling has been used to represent the groundwater reservoir in detail, simplifying the surface-water contribution only as a boundary condition of lower importance. Even in aquifers where fluxes between surface water and groundwater are limited, in the long term it is essential to consider the entire hydrologic cycle as a single system, in order to quantify the water exchanges between the different reservoirs (Mugunthan et al. 2017). This study focuses on the Firenze-Prato-Pistoia basin (Tuscany region, central Italy), which has a long experience of land subsidence and ground deformation, observed since the early 1960s (Fondelli 1975;Fancelli et al. 1980). During the last three decades, with the development of the InSAR technique, it was possible to analyze the evolution of subsidence at basin scale with millimetric precision. Canuti et al. (2006) undertook one of the first worldwide studies of this kind using remote sensing (ERS-1/2) data, showing how both the areas of Prato and Pistoia were affected by subsidence in the period 1992-2000. Prato recorded the highest subsidence rates of the area during that period, up to -20 mm/year. Rosi et al. (2016) used ENVISAT data (2001, detecting a strong increase of deformation rates, up to almost -30 mm/year, in the southeastern part of Pistoia, extending outside the historical city center towards the town of Agliana,~10 km to the south-east. The total area of subsiding 'bowl' identified had an extent of 82 km 2 ; on the other hand, as part of the same Rosi et al. (2016) study, Prato revealed a clear inversion of ground deformation, passing from a high subsidence rates to a +5 mm/year uplift. This interesting phenomenon was clearly correlated to the cessation of the local textile industry and the consequent rise of the water table and quasi-elastick rebound of the surface. Del , analyzed for the first time Sentinel-1 PSI data along the basin, detecting two main subsidence bowls close to Pistoia. The first one, located in the plant nursery area in the southeastern portion of the city, showed a mean displacement rate of -10 to -15 mm/year, slightly lower than what was detected by Rosi et al. (2016). The second subsidence bowl was detected in the historical center of Pistoia and it is characterized by an average velocity of about -15 to -20 mm/year.
The effects of subsidence on Pistoia city buildings are well known. They have been documented since the 1970-1980 decade, with the studies of Fondelli (1975) and Fancelli et al. (1980), who investigated instability phenomena recorded in Pistoia historical center during the last year of the 1960s. In recent times, Ezquerro et al. (2020) performed two damage assessment surveys in the city of Pistoia, in order to detect the extent of damage to edifices that could be potentially ascribed to the ground displacement detected by Sentinel-1 PSI in past years. The study resulted in the identification of more than 200 buildings with different levels of damage, demonstrating that the effect of subsidence in the city of Pistoia is active and cannot be neglected; however, the trigger mechanisms of the observed ground displacement were not completely investigated.
In the present study, in order to identify the causes of the subsidence phenomena that are still affecting Pistoia area, and to reduce the risk associated with their potential consequences, a numerical model of the Pistoia aquifer was developed. The model is based on an integrated simulation developed combining the US Geological Survey (USGS) finite-difference groundwater model MODFLOW (McDonald and Harbaugh 1984;Harbaugh et al. 2000) with MOBIDIC hydrologic model (Castelli et al. 2009;Campo et al. 2006). Thanks to highresolution ground displacement measurements (Sentinel-1 PSI data) and to the available hydraulic head data of the Pistoia aquifer, the groundwater model calibration was carried out, focusing on hydraulic conductivity and the elastic and inelastic skeletal storage coefficients of the aquifer. The understanding of how global changes affect water resources around the world is limited. Potential impacts of climate variation on surface water have been studied in detail by some authors (Green 2016;Shrestha et al. 2016). However, studies of how subsurface waters will respond to the combined effect of climate change and human activities have started only in the last decade (Green et al. 2011;Taylor et al. 2013;Burbey and Zhang 2015).
There is abundant evidence on the vulnerability of both surface and subsurface water resources related to climate change, including all the potential consequences for society and ecosystems (Bates et al. 2008). The assessment of climate-change effects on groundwater conditions, even if essential in a long-term hydrogeological analysis, can be very challenging, because climate variability may lead to direct and indirect repercussions for hydrogeological processes, which can be difficult to quantify (Dettinger and Earman 2007).
Increased variability in precipitation patterns and more extreme weather events caused by climate change can lead to longer periods of droughts and floods, which directly affect the availability of groundwater. In long periods of drought, there is a higher risk of groundwater resource depletion, especially in locations that depend principally on subsurface water for freshwater supply. At the same time, indirect climatechange impacts such as the intensification of human activities and/or modifications in land use, may result in the increase of groundwater demand. This can be very impactful for highly stressed aquifers, where a further increase of external stress can lead to a worsening of the aquifer condition and of related resource problems.
Land subsidence in the Pistoia area is still an active and evolving phenomenon. The analysis of potential ground displacement patterns that could affect the area in the future, as a result of the combined effect of climate variations and changing groundwater withdrawal patterns, may represent a valuable tool for urban planners and policy makers. To evaluate the effect of climate change on groundwater state, and thus on related land subsidence, a forecasting analysis using predictive global and regional climate models, in conjunction with developed numerical simulations, was carried out.
This paper represents one of the first attempts to investigate the spatial and temporal evolution of subsidence induced by groundwater-level variations, while considering the combined effects of hydrogeology, surface hydrology processes and climate-change effects at basin scale. Such effects are commonly neglected by many studies reported in the literature, but they can represent an important variable in long-term hydrogeologic simulations. Furthermore, the subsidence vulnerability maps developed by this study can represent a valuable tool for urban planners and authorities that are seeking to identify the most susceptible areas within the basin and for risk assessment analysis.

Study area
The city of Pistoia rises in the north-western corner of the Middle Valdarno basin, also known as Fi-Po-Pt basin (Firenze-Prato-Pistoia basin) in Tuscany region, central Italy. The valley is bordered by the Northern Apennines to the north and east, by the Chianti and Senese hillslopes to the south, and by the Valdinievole Valley and the Serravalle Pistoiese relief to the west. It is a 35-km wide and 100-km long intermontane sedimentary basin and has an extension of approximately 824 km 2 , with a mean elevation of ca. 50 m above mean sea level (amsl; Capecchi et al. 1976). Fi-Po-Pt basin represents one of the most important tectonic basins of the Northern Apennines, developed parallel to the main chain axis starting from the Neogene period (Boccaletti et al. 2001). The genesis of the basin is related to the extensional tectonic regime related to the opening of the Tyrrhenian Sea starting from the Upper Tortonian (Boccaletti and Guazzone 1974). The substratum of the basin is mainly composed of metamorphic Ligurian units which overlie the turbiditic formations of the Tuscan unit (Macigno sandstone). Starting from upper Pliocene, the basin is filled up with fluvial and lacustrine unconsolidated sediments with different thicknesses, reaching in the central area of the basin a depth of more than 500 m. The imposed river system strongly reflects the high anthropic impact in the area, and it is represented by the Ombrone River and its tributaries in the Pistoia province, by the Bisenzio River close to Prato, and by the Arno River in the Firenze area ( Fig. 1).
Pistoia is a medieval city with a population of about 90,000 inhabitants and an important historical center with many cultural buildings. Pistoia agricultural area (labeled "Bottegone" in Fig. 1) is characterized by the presence of several plant nurseries, representing one of the main economic activities of the area.
From a geological point of view, the area of Pistoia it is characterized by a Cretaceous-Paleogene basement constituted by a stratigraphic column starting with the sandy flysch of the Macigno Formation (Oligocene-Miocene), the Canetolo Complex clay schists (Eocene-Oligocene), the Ophiolitic Sequence (Lower Jurassic-Lower Cretaceous) and the Upper Cretaceous-middle Eocene Calvana Supergroup (Fig. 1;Capecchi and Pranzini 1986). Here, the compressible lacustrine and fluvial sediments have a thickness that increases from the reliefs surrounding the basin to more than 400 m in the southeastern part of the plain, along the axis of the basin. The alluvial deposits in the north-west of Pistoia are characterized by 30-40 m of coarse sediments, representing the fan delta of Ombrone River, where the main phreatic aquifer of Pistoia area is set. This aquifer is largely exploited for the city of Pistoia's freshwater supply. Several pumping wells are evenly spatially distributed outside of the city center, where almost no other groundwater extraction points are found, excepting for a few medieval-age dug wells. The surficial aquifer is directly connected to the surface network of Ombrone River, flowing in the western part of the city and representing its main source of recharge.
In the southeastern part of the basin, where clay and silt are more abundant, the phreatic aquifer is limited and barely exploited. Here, the main groundwater sources are represented by several multi-layers confined aquifers, set in gravelly and sandy lenses at different depths, that are characterized by high groundwater extraction rates for agricultural purposes. The sandstone rocks that underlie the basin represent the main recharge source area of Pistoia aquifer, transmitting water from surrounding reliefs through the alluvial deposit. Some additional recharge may occur from the east side of the basin, from Prato aquifer, but this amount is expected to be rather small (Capecchi and Pranzini 1986).

PSI data analysis
ENVISAT products cover the time period 2003-2010 and were processed in the framework of the PST-A (Piano Straordinario di Telerilevamento -Extraordinary Plan for Remote Sensing of the Environment) project managed by the Italian Ministry of Environment, Sea and Territory (MATTM). The PST-A project is the first worldwide example of a PSI national product, derived by analyzing ERS-1/2, ENVISAT and COSMO-SkyMed radar images in both orbits. The whole Italian territory was covered by PSI results with 28 million points for composing the final ENVISAT deformation map, with an average density of 50 PS/km 2 . The radar images were processed by means of the combined efforts of two production teams who used the PSInSAR (permanent scatterer interferometry synthetic aperture radar, Ferretti et al. 2001) and PSP (persistent scatterer pairs, Costantini et al. 2008) algorithms. Both approaches can define stable radar targets, i.e. the persistent scatterers (PS), as measurement points (MP). Each MP is characterized by a value of velocity, estimated by applying a linear deformation model, and by a time series of deformation measured along the line of sight (LoS). The density of MP mainly depends on the: (1) land cover, vegetation and wet surfaces, as well as the presence of frequent surface changes, e.g. mining activity or construction works; (2) type and style of deformation, i.e. fast motions are difficult to capture since they cause phase aliasing; and (3) topography of the area of interest, i.e. steep slopes can cause geometrical image effects such as shadowing or layover. Two high-performance computing infrastructures were used to process the ENVISAT data, each one consisting of clusters of computing nodes interconnected with large bandwidth and low latency networks (Costantini et al. 2017). The minimum processing unit was the frame, i.e. 100 × 100 km for C-band data. One single reference point for each frame was defined and selected automatically depending on the: (1) reliability of the measurement point; (2) small model residuals; (3) position within the frame; and (4) type of scatterer (e.g. a building in a historical city center was preferred over a shed in a recent urbanization). Atmospheric noise was removed using datadriven approaches based on adaptive filters (Costantini et al. 2017). The single frames were mosaicked without the support of global navigation satellite system (GNSS) data. The overlap areas between frames were used for consistency checking. The reader can refer to Costantini et al. (2017) for more detailed information regarding this wide-area processing approach.
Sentinel-1 C-band images were processed by means of a parallelized SqueeSAR algorithm (Ferretti et al. 2011). SqueeSAR is a homogeneous distributed scatterers inteferometry (HDSI) approach which relies on a new category of stable radar targets: the distributed scatterer (DS). Unlike the point-like PS, a DS is considered a distributed target whose signal is averaged from multiple targets with homogeneous signals. DSs are usually detected in uncultivated or debris-covered areas. The selection of a DS is based on the Kolmogorov-Smirnov test which is used to select statistically homogeneous families of pixels depending on the distribution of amplitude values within a predefined searching window (Kvam and Vidakovic 2007). In summary, the SqueeSAR algorithm firstly analyzes the amplitude characteristics of the single pixel of the radar image, then, through the Kolmogorov-Smirnov test, defines the families of homogeneous pixels which are going to be the DS candidate. Finally, the PS and DS are jointly processed by means of a classical PSInSAR approach. SqueeSAR allows for improvement in the accuracy of time series and velocity estimations, and it also increases the MP density in peri-urban areas and especially in mountain areas with widespread debris deposits and low vegetation cover. The data used in this work derive from a specific application of the SqueeSAR algorithm. Such data are part of the first worldwide example of a monitoring system totally based on interferometric products (Raspini et al. , 2019. The deformation maps in ascending and descending orbits derived for both ENVISAT and Sentinel-1 image stacks are the input for the estimation of the vertical component of motion which is the main calibration input for the hydrogeological model. In this work, the horizontal component of motion was not taken into account since it is negligible in the largest part of the area of interest. This is confirmed by the vertical-and horizontal-component maps included as electronic supplementary material (ESM). Only the eastern portion of the Pistoia city center records a small horizontal component in the direction of the center of the subsidence bowl, with average deformation rates below 5 mm/year. This deformation is localized, not consistent throughout the whole subsidence bowl, and it is related to lateral strain rates (Ezquerro et al. 2020). Thus, it is assumed to be not representative for the subsidence along the basin. The assumption of the vertical component as the most representative for the description of subsidence is consistent with previous PSI-related works such as Osmanoğlu et al. (2011), Higgins et al. (2014, Miller and Shirzaei (2015) and Del . The well-known methodology proposed by Notti et al. (2014) is applied to derive the north-south (vertical, V V ) component of motion. The vertical component is calculated as follows: where E and H are the direction cosines for the two orbits (a: ascending; d: descending) and vLOS is the measured velocity along the line of sight of each satellite. The cosine depends on the incidence angle (α) and on the LOS azimuth angle (γ) of the geometry of acquisition. The angles are expressed in degrees. The generic formulas for the calculation of the cosines are reported as: The calculation of the vertical component is performed using a common geographical information system (GIS). Firstly, the area of interest is subdivided in a 100 × 100 m regular grid. Then, it is verified for each cell if at least one MP for each orbit was resampled. For all the cells respecting this requirement, a synthetic MP is derived. The vertical velocity is finally computed for each synthetic MP by means of Eq. (1).
Displacement time series represent the most advanced PSI product and allow one to follow the temporal behavior of ground displacement over time. Time series are a powerful tool to verify the presence of trend changes in the normal and linear behavior of a PS point (Raspini et al. 2019). In this case, single time series were extracted in proximity of each sampling point of the model (within 50 m) and used for model validation and comparison. To avoid calculation issues due to a large amount of input data, 247 high coherence (>0.9) PS were selected, considering their position within the model domain and their reliability in terms of ground displacement time series. Time series are derived only for the MP within the search radius with temporal coherence higher than 0.9. Considering these requirements, a total of 247 MP is selected for the model calibration. If multiple MP fulfill the requirements, the average time series is calculated.

Hydrological modelling
MOdello di Bilancio Idrologico DIstribuito e Continuo, Distributed and Continuous Hydrologic Balance Model (MOBIDIC) is a fully distributed and raster-based hydrological model for water balance calculation, qualitative and quantitative water resource assessment and flood risk management (Castelli et al. 2009;Campo et al. 2006). MOBIDIC hydrological simulation was set considering, as model domain, the entire basin of Ombrone River, where the city of Pistoia rises. For the Tuscany region, quasi-continuous meteorological data time series are available from January 1992 to December 2017. The MOBIDIC structure requires that the temporal duration of each simulation is defined by the hydrometeorological data (precipitation, maximum and minimum temperature, wind speed, solar radiation and air humidity) time series length: 26 years in this case . The first 2 years of simulation (1992)(1993) were used to tune the model, in order to give enough time to the model to fill up its reservoirs and to reach stable initial values of hydrologic variables. The calibration phase involved 5 years of the simulated temporal domain (2010)(2011)(2012)(2013)(2014)(2015), chosen for its good availability of continuous field data, in a period as close as possible to the maximum magnitude of the observed subsidence events.
Calibration of the MOBIDIC global hydrologic semi empirical parameters was carried out by means of the calibration module integrated in the code. This module searches the minima of the objective function given by the weighted sum, for each available hydrometric series, of the normalized differences between measured and computed data (Eq. 4). In the objective function, four different parameters can be considered: discharge, cumulative flow volume, flow duration and peak discharge, and for each term the relative weight in the objective function must be specified: where n is the number of parameters involved in the calibration phase, R i is the residual error for each parameter and w i is the weighting factor attributed to each parameter. For the calibration phase, two measurement stations (Pontelungo and Poggio a Caiano, Fig. 1) located within the Ombrone watershed and quasi-continuous-recorded river discharge data for the whole calibration period (2010)(2011)(2012)(2013)(2014)(2015) were considered. The focus of the calibration phase was to identify the optimal values for the four input hydrologic parameters of MOBIDIC, in order to reach the best fit between model values and observations. The initial values for the global hydrologic calibration parameters were obtained from the Arno River hydrological model, carried out with MOBIDIC by Castelli et al. (2009), and currently in use as a flood risk assessment tool by the Arno River Basin Authority Hydrogeological modelling MODFLOW is the USGS three-dimensional (3D) finitedifference modular groundwater model (McDonald and Harbaugh 1984;Harbaugh et al. 2000). In MODFLOW the groundwater flow equation is solved using a block-centered finite-difference approach, subdividing the horizontal modelled domain into a grid constituted of an arbitrary number of rectangular cells and the vertical domain divided into different layers, representing the geological variability.
The horizontal discretization of Pistoia aquifer was modelled by simplifying the real geometry of the aquifer, characterized by a high level of detail that could had led to convergence issues in the MODFLOW simulations (Fig. 2a). The model area is 18.0 × 18.0 km 2 and it was discretized into a grid of 360 × 360 cells using a uniform cell size of 50 m.
The vertical domain was discretized into three model layers that were used to represent three aquifer zones (labelled 1, 2 and 3 in b). Model layer 1 corresponds to the near-surface phreatic aquifer, directly connected with the Ombrone River network, and characterized by an average thickness of 8-10 m all over the plain. An exception is represented by the layers of thick gravel and pebbles in the northern portion of the area, coinciding with the Ombrone River fan delta. Beneath the surficial aquifer, alluvial deposits can be divided into two sublayers with different sedimentary facies (Capecchi et al. 1976;Capecchi and Pranzini 1986). The deeper layer is mainly constituted by upper Pliocene lacustrine sediments, generating sand lenses embedded into a finer matrix that limits the aquifer permeability. The shallower layer is comprised of coarser Pleistocene materials, where the fluvial component becomes predominant and constitutes a layer of higher permeability. This part of the vertical model domain is simplified by means of two layers, representing the middlezone (layer 2) and the deep-zone (layer 3) aquifers. The bottom of layer 2 was estimated based on borehole data available for the study area, while the bottom of layer 3 was set in correspondence with the stratigraphic contact between the basement and the valley fill.
The model consists of 25 stress periods, composed by a stationary stress period and 24 transient ones, corresponding to the 1994-2017 period, with year as a unit time step. The lengths of the stress periods varied from 10 days for the stationary stress period to 365 days for the transient stress periods.
The main water input for the Pistoia phreatic aquifer is the surficial network, providing a constant source of water especially during the rainiest months of the year. The river system (Fig. 2a) was modelled by means of the river boundary conditions (RIV package), obtaining the water stage values from the MOBIDIC hydrologic simulation.
Recharge to the confined unit is quite limited in short-term periods and can be identified in the percolation water that flows from the high-permeability sandstone formations that underlie the basin into the aquifers. This contribution was taken into account in layers 2 and 3 by setting injection well cells (modelled with MODFLOW WEL package) all around the model domain (Fig. 2a). Another low-entity recharge source for the surficial layers is represented by the direct precipitation on the valley floor, modelled with the RCH (recharge) package. Both contributions were directly quantified thanks to the percolation variable simulated with MOBIDIC code.
Some additional water exchanges, in both directions, may occur in the east side of the basin, where a real physical limit with Prato aquifer does not exist (Fig. 2a). This small amount of water was considered by setting a general head boundary condition (GHB package) along the east side of the model domain and its value was set between 40 and 80 m amsl, basing on available hydraulic head observations.
Outflows from the aquifer system are mainly by groundwater extraction for drinking water supply and irrigation purposes, especially for the numerous horticulture activities (flower and plant nurseries) present in the south-east area of Pistoia. Groundwater extraction occurs from a total of 5,910 extraction points (2,397 for the surficial aquifer, 2,844 for middle-zone aquifer and 669 for deep-zone aquifer) with a total pumping rate of about 50 Hm 3 /year. Detailed information about groundwater extraction (well depth, location and discharge) were obtained from the online database of the local authorities (Arno River Basin Authority 2013). Wells were modelled by means of the WEL (well) package.
Since the basin sediment fill is interrupted by the surrounding less-permeable reliefs, the lateral boundaries of the model are all modelled as no-flow conditions in all three layers.
For the simulation of ground displacement related to groundwater withdrawal, the model uses the Subsidence and Aquifer-System Compaction Package (SUB package) included in the MODFLOW code (Höffmann et al. 2003). Key parameters considered by MODFLOW that affect both hydraulic head and subsidence distribution all over the aquifer are hydraulic conductivity (K) and elastic (S ke ) and inelastic (S kv ) skeletal storage coefficients. Such parameters, together with groundwater extraction rate, are largely responsible for the observed land-surface deformation and water level distributions, but the assessment of their spatial variability and distribution can be very challenging and must be properly calibrated during the development of the model (Yan and Burbey 2008).
Hydraulic head maps of the Ombrone aquifer that referred to 1994 and 1996 were used for steady state calibration of the model; these periods represent a quite stationary period with respect to groundwater level, with very low variations registered in the measured hydraulic head (Fig. 2a).
Piezometric levels are the most popular type of observation data used for parameter estimation in groundwater flow models (Poeter and Hill 1997). However, in regions where a complete network of hydraulic head monitoring does not exist, obtaining enough data can be problematic for economic reasons. Satellite interferometric data support the measurement of land subsidence at basin scale and can be used together with sparse drawdown or hydraulic head data to improve the estimation of calibration parameters (Yan and Burbey 2008, Bell et al. 2008, Burbey and Zhang 2015. Since no continuous hydraulic head data were available for the present study area to calibrate the transitory groundwater model, highdetail Sentinel-1 ground displacement data for the 2015-2017 period were used for the purpose. Additionally, ENVISAT PSI data acquired between 2003 and 2010 were compared with simulation results, in order to validate the numerical model and to confirm its reliability during the whole analyzed temporal domain. The preconditioned conjugate-gradient (PCG) method was used as the numerical solver, with HCLOSE (head change criterion for convergence) and RCLOSE (residual criterion for convergence) values set to 0.001 and 1, respectively. The numbers of outer and inner iterations were set to 200 and 30, respectively.

Future climate scenarios
Global climate models (GCMs) and regional climate models (RCMs) are some of the most powerful tools for simulating future climate scenarios, providing all the basic information needed for climate-change impact assessment at many different scales (Fang et al. 2015). The fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change released a new set of scenarios for climate-change forecasting, introducing the representative concentration pathways (RCPs; AR5 2014). Four RCPs were introduced, called respectively RCP2.6, RCP4.5, RCP6 and RCP8.5, which get their names from the possible values of radiative forcing that they are simulating for the year 2100. The coarse resolution of GCMs makes them unable to represent climate variations at regional or finer scale (Ramirez-Villegas et al. 2013). Regional climate models guarantee higher resolution than global climate models, preserving the coherence between the meteorological variables' distribution and the topographic characteristics of the Earth's surface. However, RCMs are still affected by many systematic errors (bias) caused by their rough resolution and by the assumptions made during their development (Christensen et al. 2008;Suklitsch et al. 2011;Kjellström et al. 2010). In order to perform reliable simulations, the raw data from GCMs and RCMs need to be biascorrected before they are used as input for a climate forecasting model. Many different bias correction procedures have been proposed in the literature to transpose climate model output to local scale, from simple downscale methods to more sophisticated geostatistical approaches (Teutschbein and Seibert 2012). In the present study, the linear-scale method (LS) was used (Ines and Hansen 2006). This method is based on a constant correction factor given by the differences between GCMs/RCMs output and observations for each month in the historical reference period (Eq. 4). The correction factor is then applied to the raw climate model data in an additive or multiplicative way, in order to obtain bias-corrected time series of the studied variables (Eq. 5).
where T cor,m,d represents the bias-corrected RCM data in the day d of the month m, T raw,m,d represents the raw RCM data in the day d of the month m, μ(T obs,m ) represents the mean value of observations at the month m in the reference period, and μ(T ref,m ) represents the mean value of observations at the month m in the reference period. The LS method is capable of perfectly matching the monthly mean of the corrected RCM output with that of the observed values (Lenderink et al. 2007).
The corrected outputs of the global and regional climate models were used as input parameters for the developed hydrologic model, in order to forecast the response of the hydrologic and hydrogeologic states of Ombrone basin to the climate variation that is expected to occur up to 2050.

PSI velocity maps
The velocity maps obtained for the ENVISAT and Sentinel-1 datasets in descending orbit are reported in Fig. 3a, b. The velocity maps for the ascending orbit are shown in figures S1 and S2 of the ESM. The stability threshold (±2 mm/year) is defined as two times the standard deviation of the velocity values.
The ENVISAT velocity map (Fig. 3a) shows a clear subsidence area extending between the southern part of Pistoia, outside of the historical city center, to Agliana town. The highest deformation rates are recorded in the area of Bottegone, which is the center of the plant nursery activity area and where the groundwater exploitation is more intense. In this area, LOS velocities reach values up to 30 mm/year. A localized subsidence bowl is located near Agliana municipality, this time probably related to the industrial use of groundwater.
The Sentinel-1 velocity map (Fig. 3b) presents some evident variations: (1) the main subsidence bowl elongated along the basin axis is reduced in extent with Bottegone as its center, with LOS velocities comparable to the ENVISAT period; (2) subsidence rates around the city of Agliana town are greatly reduced with some sectors below the stability range and (3) a new subsidence bowl located in the city center of Pistoia with LOS velocities up to -20 mm/year. The reduction of the extension of the main subsidence bowl is related to the reduced activity of the textile industry in the area of Prato-Agliana (Rosi et al. 2016) and to a better use of the underground water resources in the plant nursery area. Figure 3c, d presents the vertical velocity (V V ) maps for the two analyzed datasets. The ENVISAT and Sentinel-1 V V maps confirm the deformation patterns highlighted by the LOS velocity maps. Figure S3 of the ESM presents the time series of deformation for points 1 and 2 in Fig. 3, referring to the area of Bottegone and Pistoia city center, respectively. Time series were produced for both orbits and include both interferometric datasets in order to highlight any trend change. The time series for Bottegone show a linear trend over the whole period with a more evident 6-month seasonality in the ENVISAT period. A reduction of the accumulated deformation is evident in the Sentinel-1 period. The time series of Pistoia well explain the occurrence of the new subsidence bowl. The acceleration between the ENVISAT and the Sentinel-1 period is clear.
The main characteristics and some quality parameters of the interferometric datasets used in this work are presented in Table S1 in the ESM. It is worth nothing the increase of image availability of Sentinel-1 with respect to the older generation C-band satellite. This great data availability permits one to increase the quality of satellite-derived estimation and grants an update frequency which is optimal to follow a phenomenon such as subsidence related to groundwater withdrawal. The standard deviation of the velocity values is below 1 mm/year, testifying the quality of the interferometric data.

Hydrological modelling
Model calibration was performed using the MOBIDIC calibration module, attributing a different weight to the four parameters considered by the objective function. To enhance the importance of flow duration over long periods of time, an 85% weight was given to the flow duration parameter, while a 5% weight was attributed to discharge, cumulative flow volume and peak discharge. Calibrated global hydrologic parameters and their relative initial values are presented in Table 1: To evaluate the goodness of the model calibration, duration curves of Ombrone River at Poggio a Caiano gauging stations are created (Fig. 4). The resulting hydrologic model is very reliable, with flow duration curves that match with the observations in all the years of simulation. During the calibration phase, the 85% weight in the objective function was attributed to the flow duration parameter, in trying to reduce the difference between observed and modelled duration curves as much as possible. Only a 5% weight was given to discharge, cumulative flow volume and flow peaks, resulting in major The analytical evaluation of the model reliability was carried out by means of the Nash-Sutcliffe efficiency index (NSE), the root-mean-square error (RMSE) and the normalized root-mean-square error (NRMSE) statistical indexes (

Hydrogeological modelling
Groundwater model calibration was carried out by means of the PEST code (Model Independent Parameter ESTimation) which was used together with MODFLOW-2000 to conduct an inverse modeling procedure (Doherty 2010). Since only scarce data and/or punctual data on the main hydrodynamic parameters were available in the area, the main goal of the calibration phase was to estimate hydraulic conductivity and  Hillslope flow coefficient 5.00 × 10 −5 7.77 × 10 −5 elastic and inelastic storage coefficients for the Pistoia multilayered aquifer. The starting values of such parameters were considered to be constant for each model layer and were obtained from the analysis of the few field data available on the databases of local authorities (Table 3). Thanks to the availability of high-resolution PSI data, model calibration allowed for definition of the spatial distribution of the hydrodynamics parameters within Pistoia aquifer, in order to take into consideration the influence of their distribution in controlling subsidence phenomena in the study area. The anisotropy of the hydraulic conductivities is specified as K h /K v =10, where K h and K v are the horizontal and vertical hydraulic conductivities, respectively.
Since the interbeds included in layers 2 and 3 of Pistoia aquifer can be several meters thick, the assumption that heads within interbeds equilibrate instantaneously when a change in aquifer head occurs cannot be made, and delay interbeds must be considered for subsidence modelling (the equilibration of hydraulic heads within the interbeds is delayed if compared with head changes in the aquifer). On the other hand, interbeds contained in layer 1 are less thick (decimetric scale) all over the model domain. For such a case, no delay interbeds are adopted. Simulated water levels and ground displacement values were compared with available hydraulic head data and PS subsidence data all over Pistoia aquifer, in order to evaluate the accuracy of the groundwater model calibration (Fig. 5). The water balance of the first stationary time step is also provided. Residual mean, residual standard deviation, root means square error (RMS) and determination coefficient R 2 were used to evaluate the model fit. Figure 5b shows the good fit (R 2 = 0.987) obtained between modelled and observed hydraulic head values after the calibration phase, considering their comparison with precalibration parameter (R 2 = 0.745, Fig. 5a). A total of 215 points were used to derive the graph. Model residuals, calculated as the difference between model values and observations, are on average equal to -0.23 m, the average standard deviation is 1.38 m and the scaled RMS is 2.5%. Assuming a piezometric variability of 56 m in the hydraulic head observations, both errors are below 10% of the total range of observations; this is an acceptable error for a long-period regional model.
As for the hydraulic head, the comparison between simulated and measured ground deformation data testifies to the goodness of the calibration, with a R 2 value that rises from 0.640 (considering starting values of calibrated hydrodynamics parameters) to a value of R 2 near 1 (0.9784, Fig. 5c, d). Residuals show an average value of 0.6 mm, with a standard deviation of ±5 mm. Considering a difference of about 50 mm between the maximum and minimum observed subsidence, the model residual falls within the 10% of acceptable error, confirming the reliability of the model.
The water balance of Pistoia aquifer during the stationary stress period is shown in Fig. 5e. Aquifer recharge, occurring as the direct effect of precipitation, and river inflow and outflow are associated only with the surficial phreatic aquifer. On the other hand, groundwater fluxes from surrounding aquifers (simulated as GHB boundary condition) and lateral affluxes represent the main inflows of the deep aquifers.
The final calibrated values of elastic and inelastic storage coefficients for layer 2 of the numerical groundwater model are provided in Fig. 5f. Higher calibrated values of S ke and S kv characterize the area of the city of Pistoia and of Bottegone, suggesting that fine compressible interbeds may be abundant in these locations, representing one of the main factors that control subsidence in the area. Very similar storage coefficient distributions characterize layers 1 and 3.
Time series of simulated land subsidence in the Pistoia historical center and plant nursery area are presented in Fig.  6, together with the vertical component of ground displacement obtained from Sentinel-1 PSI data for the 2015-2017 calibration period.
By means of the millimetric accuracy and the weekly acquisition frequency, Sentinel-1 time series describe a Table 3 Starting values of groundwater model calibration parameters: horizontal hydraulic conductivity (K x , K y ), specific storage (S s ), porosity (n) and elastic and inelastic skeletal specific storage coefficients (S ke and S kv ) 1 2 . 3 × 1 0 −3 2.3 × 10 −4 0.72 0.25 1.0 × 10 −6 1.0 × 10 −4 2 6 . 1 4 × 1 0 −5 6.14 × 10 −6 0.02 0.25 3.0 × 10 −6 3.0 × 10 −4 3 1 0 −6 10 −7 0.01 0.25 6.0 × 10 −6 6.0 × 10 −4 Nash-Sutcliffe efficiency index (NSE), the root-mean-square error (RMSE) and the normalized root-mean-square error (NRMSE) clear seasonal trend within a 6-month period, except for the time series shown in Fig 6b, which registered a bilinear trend with high periodic fluctuations (intrinsic for the PSI technique). The time series recorded an average velocity in the 3-year-long monitored period between -13 (Fig. 6c) and -16 mm/year (Fig. 6d) and a total displacement between -40 and -50 mm. For those time series showing a clear 6-month variation, subsidence is registered during the dry season and a lower-magnitude recovery during the wet season. The maximum subsidence reaches values of -25 mm in the dry season, whereas the maximum recovery is equal to +10 mm in autumn and winter. Because of the yearly resolution of the groundwater model, time series of simulated subsidence do not show any cyclical short-term oscillations but confirm the general trend of deformation, which is comparable with observed data values. In 2016, for both plant nursery areas and for the western part of the historical center, simulated ground displacement shows a stabilization in the general deformation pattern (Fig. 6a-c). Such differences with observations are related to a temporary recovery of the hydraulic heads simulated by the groundwater At the end of the simulated period, the model is able to reproduce the magnitude of observed subsidence very well, showing a good agreement with PSI data. The comparison between modelled and observed cumulated subsidence confirms the goodness of the calibration phase, with only small differences, lower than 10%, observed in the total ground displacement. Only in the eastern part of the plant nursery area, the modelled total displacement exceeded the 10% error, reaching a value of 13.5% (Fig. 6c).
The subsidence model validation was carried out, for the 2003-2010 period, by means of the available ENVISAT PS data available all over Pistoia aquifer. Time series of observed and modelled ground displacement for the validation period of both the Pistoia historical center and plant nursery area are presented in Fig. 7.
ENVISAT PSI time series referring to the plant nursery area show a constant and linear trend, characterized by a mean velocity higher than -25 mm/year and a maximum cumulated displacement of about -240 to -250 mm. Despite the monthly or bimonthly acquisition time, ground displacement observed by ENVISAT sensors is characterized by a clear 6-month seasonal oscillation which overlies the general linear trend. The subsidence model is not able to describe such a clear trend, but it can reproduce the general pattern and magnitude of ground displacement in quite good agreement with observations. For the Pistoia historical center, the ENVISAT data detected a mean velocity lower than -5 mm/year and a total displacement of -40 to -45 mm at the end of the monitored period. There is a correspondence between the model and PSI data for the western portion of the historic city center, whereas the model overestimates the total displacement in the eastern portion of the city center of about -15 mm. It is worth noting that the time series referred to this sector of the city is noisy, with strong variations probably not due to ground motion.

Subsidence forecasting
Calibrated RCM output for precipitation, maximum and minimum temperature, wind speed, solar radiation and air humidity were used as input data for the calibrated MOBIDIC model, in order to extend the hydrologic simulation of the Ombrone River network up to 2050. To certify the goodness of hydrologic simulations performed with the RCM output, the mean monthly discharge values of all RCM models were compared with the original calibrated model for the 1994-  (Fig. 8).
According to the bias correction method that was used to calibrate RCM meteorological data, a good match between mean monthly values is observed. RCM mean monthly calculated discharges are in very good agreement with calibrated MOBIDIC model output, showing similar values for each time step. Both Poggio a Caiano and Pontelungo stations show high coefficient of determination values of 0.98 and 0.96, respectively, certifying the reliability of RCM simulations.
In order to identify the potential evolution of landsubsidence patterns affecting Pistoia area in the future, a groundwater model composed by 58 stress periods was set up, including a 25-year model validation period (1992-2017) and a 33-year forecast period (2019-2050). Groundwater simulations are based on the calibrated MODFLOW model of Pistoia aquifer, using, for the definition of boundary conditions, calculated river discharge and percolation variables obtained from RCM hydrologic simulations. Three different pumping scenarios were investigated considering: (1) the current pumping rates and water dynamics, (2)  an average increase in pumping rate of +1%/year, starting from 2020 (+10% at the end of 2030, +20% at the end of 2040 and +30% at the end of 2050) and (3) an average decrease in pumping rate of -1%/year, starting from 2020 (-10% at the end of 2030, -20% at the end of 2040 and -30% at the end of 2050).
To evaluate the goodness of the groundwater simulation developed with RCM boundary conditions, a comparison between the modelled ground displacement rate and the Sentinel-1 PSI data was carried out, considering the average subsidence rate identified during the 2015-2017 period (Fig.  9). The RCM groundwater model reproduces the general shape of the subsidence bowls with quite good accuracy, both in Pistoia historical center and in the plant nursery area. Here, the model slightly underestimates the observations value, identifying a maximum ground displacement rate of about -15.5 mm/year (against 20 mm/year of PS data) and a subsidence bowl area with less spatial extent. Major discrepancies are encountered east of the plant nursery area, where the RCM model identified no subsidence occurring in the 2015-2017 period (lower than 2.5 mm/year), whereas the Sentinel-1 PSI data detected subsidence rates up to -15 mm/year.
The receiver operating characteristic (ROC) curve of the forecasting subsidence model follows the left side and the top border of the chart, demonstrating the high reliability of the predictive model. The value of the area underlying the curve (AUC) ranges between 1 and 0.5 (represented by the dotted line in Fig. 9b) and it provides a further indicator of the consistency of the predictor. Following the performance rating proposed by Swets (1988), the resulting model was mildly Fig. 9 a Pistoia aquifer ground-displacement velocity comparison between Sentinel-1 PS-InSAR data (point features) and RCM groundwater model estimations (isolines) for the 2015-2017 period, and b the receiver operating characteristic curve of the forecasting subsidence model. Forecasting model performances were analytically assessed by means of a receiver operating characteristic curve, also known as a ROC curve (b). The ROC curve represents a powerful graphical plot that illustrates the diagnostic ability of multiclass classification problems by plotting true positive rate against false positive rate of the prediction output accurate, with an AUC value of about 0.86. The analysis of the potential ground displacement velocity that could affect Pistoia aquifer in the future was performed in a decennial time frame, considering the average subsidence rate of the 2020-2030 (I), 2031-2040 (II) and 2041-2050 (III) periods (Table 4).
Assuming extension of the current water extraction rate from the aquifer, during the 2020-2030 period the evolution of subsidence in Pistoia area shows different trends depending on the location within the basin (Fig. 10a). In the Pistoia historical center, the ground deformation is around -5 to -10 mm/year, exhibiting a decrease from the measured subsidence rate at the end of 2017 (-15 to -20 mm/year). The plant nursery area keeps showing an average subsidence rate of approximatively -15 mm/year, with peaks of -20 to -30 mm/year. Forecasts for 2030-2040 indicate a slight increase in subsidence rate, with an extension of the involved area in the southeastern part of the city (Fig. 10b). Finally, at the end of the 2040-2050 period, the subsidence rate is expected to slow down in the historical center of Pistoia, while larger ground displacement is still affecting the plant nursery area, in accordance with the subsidence trend observed during the previous decade (Fig. 10c).
Considering an increase of +1%/year in pumping rate, starting from January 2020 (for an overall +10% at the end of 2030), the model results indicate that the current subsidence bowls in Pistoia city and its surroundings will expand, and the land surface of these areas will experience greater subsidence rates, especially in the agricultural areas. According to the forecasting model results, the ground displacement rate is expected to rise up to -25 mm/year in the plant nursery area and up to -15 mm/year in downtown Pistoia. During the 2030-2040 period, subsidence bowls keep expanding in all directions, involving the northwestern part of the city of Pistoia and the southeastern part of the plant nursery area. The land surface keeps subsiding rapidly, with greater velocity than in the previous decade. At the end of 2050, all investigated areas are subsiding rapidly, assuming very high velocities in both the historical center and plant nursery area. Except for a small area north-east of Pistoia, the area affected by subsidence phenomena remains constant, if compared to the previous decade. Results of subsidence forecasting analysis for the increasing pumping rate scenario are displayed in Figure S1 of the ESM.
Model output of the decreasing pumping rate scenario for 2020-2030 shows a reduction in the spatial extent and intensity of subsidence all over the study area, with velocities fixed below -5 mm/year in most parts of the region, and only a few areas are characterized by higher velocities (up to -10 mm/ year, Figure S2 in the ESM). At the end of the 2030-2040 period, the results indicate a total absence (lower than 2.5 mm/ year) of ground displacement related to groundwater exploitation for the entire basin. Only small subsidence bowls are still visible in the central part of plant nursery area, but they are limited in terms of spreading and intensity. This represents a clear indicator of the influence of groundwater-level variation in controlling aquifer compaction in alluvial environments characterized by the presence of fine compressible sediments; this emphasizes the importance of correct management of the groundwater resource. A further decrease in groundwater extraction rate, as simulated in the 2040-2050 period, does not strongly affect the residual subsidence pattern still occurring in the plant nursery area. Here, subsidence bowls are very limited in extent and they are characterized by velocity values lower than 2.5 mm/year. This is probably linked to the stabilization of water levels within the aquifer that started during the previous analyzed decade, as a direct consequence of the reduction of water extraction from the groundwater system.
Due to the absence of continuous and detailed hydraulic head data in Pistoia aquifer in recent years, a direct comparison between ground displacement rate and water-level changes cannot be performed in the present study. Updated hydraulic head measurements would be very useful to confirm the correlation between subsidence and groundwater withdrawal, also providing important information for authorities and policy makers to improve understanding of such an impactful phenomenon.  I  II  III  I  II  III  I  II  III Historical center (

Conclusions
Fi-Po-Pt (Firenze-Prato-Pistoia) basin experienced a long history of subsidence phenomena, clearly identified since the early 1990s by means of the first available InSAR data. In the present work, a hydrogeologic model of Pistoia aquifer was developed by means of the MOBIDIC-MODFLOW integrated procedure, in order to characterized subsidence patterns induced by groundwater withdrawal and aquifer compaction. Model calibration and validation were carried out thanks to high-resolution subsidence observations derived from PSI data, focusing on hydraulic conductivity and skeletal storage coefficients of the aquifer, demonstrating the benefits of using detailed subsidence data in the characterization of aquifer properties when detailed and complete hydraulic head measurements are not available. Calibrated values of storage coefficients may suggest that the combined effect of groundwater extraction and the presence of compressible fine sediments may represent the driving force of subsidence detected in Pistoia area. An evaluation of potential anomalous subsidence patterns that could be applicable to the study area in the near future was performed, using forecasted meteorological data of global and regional climate models as input for the MOBIDIC-MODFLOW modelling procedure. Ground displacement simulations were extended up to 2050, considering three different pumping-rate scenarios. This led to the development of several subsidence hazard maps of the city of Pistoia, which may represent a valuable tool for urban planners to identify the most susceptible areas to subsidence and to develop new strategies to reduce the associated risks.
Since aquifer compaction caused by groundwater withdrawals depends on several variables and factors, MODFLOW code must perform with some assumptions and limitations with respect to the discretization of the natural system. In the MODFLOW SUB package, the skeletal specific storage values, S ske and S skv , are stress dependent. As head declines and effective stress increases, these values should become smaller. As a direct consequence, calculated compaction and storage change are smaller for equivalent changes in effective stress. Unfortunately, the SUB package does not consider reductions in skeletal specific storage for subsequent stress periods and projections of future subsidence might be overestimated. Moreover, since geostatic load variations are not considered by the SUB package, only changes in effective stress caused by changes in fluid pore pressure can be simulated. Geostatic load changes may occur for the presence of new buildings or engineered structures and/or by changes in water levels in overlying unconfined aquifers, but such variations are not considered in groundwater simulations. Since neither interbeds nor model-layer thickness are adjusted to account for their compaction during the simulation, the SUB package also assumes that compaction of individual interbeds is small if compared with their total original thickness. Layer thickness is fixed during the Climate variation effects do not seem to have a strong influence on subsidence control, as ground displacement rate remains quite stable during the projected 2020-2050 period, if one considers the extension of current pumping rates in the future. Following the expectations, maximum subsidence values are forecast using the increasing pumping rate scenario, while lower subsidence values are obtained by decreasing the water extraction rate. The present work confirms the role of groundwater extraction as one of the primary driving force of subsidence patterns in the Pistoia area, showing a higher influence than short-term climate changes on controlling ground displacement.
In addition, forecasting analysis showed that if pumping from the aquifer continues uncontrolled, anomalous subsidence rates affecting the city of Pistoia are expected to persist, leading to severe consequences for buildings, engineering structures and human activities. As a first step, given the high compressibility nature of the alluvial sediments in the study area, subsidence mitigation actions should include a correct and sustainable use of the groundwater resource. New and updated water level measurements could provide useful information to assess current groundwater states in Pistoia aquifer, allowing the development of a finer and better-calibrated numerical model with improved reliability. Considering the anomalous ground-surface displacement velocity observed by PSI in the last 30 years in Pistoia surroundings related to high withdrawals rate, an adapted and forward-looking groundwater management scheme in the area is of vital importance in order to guarantee groundwater quality and quantity, and its renewability.