Impacts of 25 years of groundwater extraction on subsidence in the Mekong delta, Vietnam

Many major river deltas in the world are subsiding and consequently become increasingly vulnerable to flooding and storm surges, salinization and permanent inundation. For the Mekong Delta, annual subsidence rates up to several centimetres have been reported. Excessive groundwater extraction is suggested as the main driver. As groundwater levels drop, subsidence is induced through aquifer compaction. Over the past 25 years, groundwater exploitation has increased dramatically, transforming the delta from an almost undisturbed hydrogeological state to a situation with increasing aquifer depletion. Yet the exact contribution of groundwater exploitation to subsidence in the Mekong delta has remained unknown. In this study we deployed a delta-wide modelling approach, comprising a 3D hydrogeological model with an integrated subsidence module. This provides a quantitative spatially-explicit assessment of groundwater extraction-induced subsidence for the entire Mekong delta since the start of widespread overexploitation of the groundwater reserves. We find that subsidence related to groundwater extraction has gradually increased in the past decades with highest sinking rates at present. During the past 25 years, the delta sank on average ∼18 cm as a consequence of groundwater withdrawal. Current average subsidence rates due to groundwater extraction in our best estimate model amount to 1.1 cm yr−1, with areas subsiding over 2.5 cm yr−1, outpacing global sea level rise almost by an order of magnitude. Given the increasing trends in groundwater demand in the delta, the current rates are likely to increase in the near future.

The following document provides additional description and presentation of data and methods as referred to in the manuscript.

S1. Subsurface model
We adopt the geological subsurface schematization of the Mekong delta (MKD) and the Saigon delta made by the Division for Geological Mapping for the South of Vietnam (DGMS, 2004). As these two delta systems share the same depositional basin, their deposits and groundwater systems are interconnected. Ten hydrogeological cross-sections based on 95 linearly interpolated, geological borehole logs were used to construct the 3D subsurface model of MKD using a linear interpolation of iMOD SolidTool (Vermeulen et al 2016). The deepest confined aquifer (Upper-Middle Miocene age) was not represented in the model due to data paucity. This aquifer is less relevant for this study, as no groundwater extraction is reported for it. Surface elevation was based on a detailed digital elevation model of the Mekong delta supplemented with Shuttle Radar Topography Mission (SRTM) elevation data for areas outside the delta. This resulted in a 15 layered, subsurface model, representing 7 aquifers, 7 aquitards and a phreatic top layer (Table S1). The minimum layer thickness was set to two meters, as required for numerical functioning of the current build of the subsidence module in iMOD (SUB-CR). Furthermore, an adjustment was made in Soc Trang province: an incorrectly interpreted aquitard thickness in a geological borehole log was manually corrected to a minimum thickness of 5 meters. In the current subsurface model, ~60% of the total volume is represented as aquitard (clay and silt) and ~40% as aquifer (sand). A lithological analysis of over 700 core logs revealed a sand versus silt and clay presence for the Mekong delta of respectively 60% and 40%. This implies that about 1/3 of the aquitards in the current subsurface model actually consist of sand, which is reflected by the abundant sand lenses present in the core logs. This is taken into account in the subsidence calculations to avoid overestimation of subsidence. Table S1. Geological and hydrogeological units of the MKD and corresponding model discretization.
In total ~40% of the subsurface presented in the model is aquifer and ~60% is aquitard. Model layer 1 is not shown in this table, as it is not based on a specific geological formation, but represents the phreatic top layer at the delta surface.

S3.1 Mekong delta
For the MKD two data sets of past groundwater are used, reporting well location (x-and ycoordinates), exploited aquifer and daily extracted volume (DWRPIS, 2010). The first dataset reports on large, industrial extraction wells with an extraction permit (>200m 3 /day) covering the start of the modelling period until 2011. The year of registration is assumed to be the starting point of the extraction, and the permitted volume to be the actual daily extraction value. The second dataset is based on a large delta-wide survey held in 2010 using household density and interviews to estimate well depth and daily extraction by the Division of Water Resources Planning and Investigation of South Vietnam (DWRPIS). The outcomes of the survey were grouped in fictive well locations randomly distributed over freshwater zones for each province with extraction rates to match total reported provincial extractions for each aquifer (Table S3). We assumed that the unregistered household groundwater equally grew with the industrial demand. Therefore, the growth recorded in the registered industrial extractions was used to extrapolate the 2010 survey results back to 1991.
Following 2011, an annual increase of 2.5% was assumed. This results in a total extraction in the MKD approaching 2.5 million m 3 /day at the end of the modelling period (Fig. 3, main article). The Lower Pleistocene (48% of the total volume) and the Lower Pliocene (26% of the total volume) aquifers are the most heavily-exploited in the MKD ( Figure S1).

S3.2 Ho Chi Minh province
Three datasets were combined to determine the extraction of groundwater in the HCMC province.
The first dataset reports small extractions (<10 m 3 /day) per aquifer but lacked location information.
The second dataset documented larger extractions (>10 m 3 /day) per aquifer including well coordinates. Subsequently, assuming similar extraction patterns, these well locations were adopted for the first dataset (  (Table S3).

S3.3 Uncertainty in extraction data
We identify four caused of uncertainty in the developed extraction database: 1) The spatial coverage of well data is not all-inclusive (Fig. S2). Several areas where groundwater extraction is suspected are A potential solution for the above-mentioned data deficits is to include extracted volume in the model calibration process, link extraction volume to land use practice and well locations to population density and growth over time. However, adjusting the dataset in this manner would introduce different errors and uncertainties, and requires a lot of detailed local data, as, for example, the main pumping stations producing fresh water for the city of Tra Vinh are located many kilometres away from the city.

S4.1 Initial parameterization prior to calibration
Initial values of horizontal hydraulic conductivity (K h ) for the aquifers range from 8.0 to 22.8 m/day and were derived from 999 pumping tests throughout the Mekong and the Saigon delta (DWRPIS, 2010). Aquitard layers were initially parameterized with a K h value of 0.001 m/day. No measured data was available on the vertical anisotropy of hydraulic conductivity (K h /K v ). We assumed a vertical anisotropy of 3, as general values of vertical anisotropy for aquifers range between the 1 and 3, with outliers up to 8 (Carlson, 2000). Initial specific storage coefficient (SS c in m -1 ) values for individual aquifers (ranging from 1.4 x10 -3 to 5.6 x10 -3 m -1 ) and aquitards (ranging from 3.8x10 -5 to 2.2 x10 -3 m -1 ) are taken from reported values by Haskoning et al (1999) and Giao et al (2015).

S4.2 Parameter optimisation
The automated parameter estimation in iMOD (i.e. iPEST, see Vermeulen et al 2016 for a full description of the method) was applied using piezometric head measurements from 91 monitoring wells using ten pilot points (Fig. S3). At every pilot point, K h and SS c were systematically adjusted during the PEST procedure (when total model sensitivity to the parameter > 0.5%) and interpolated (simple kriging) over the entire model extend. This was done for each model layer. The maximal multiplication factor allowed for initial values for K h and SS c was set to 200 and 1000 for respectively the aquifer and aquitard layers to ensure realistic values. The parameter set creating the smallest error between modelled hydraulic head with measured head time series between 1995-2015 was determined through numerous iterations with each measured time series equally weighted (see also  (Table S4).

Figure S3
Locations of observation wells and pilot points used during the automated parameter estimation (iPEST).

Table S4
Calibrated values for horizontal hydraulic conductivity and specific storage coefficient following the automated parameter estimation.

S5.1Determining the geotechnical parameters
The geotechnical parameterization for the abc model of the MKD subsurface was based on general relationships existing among compression parameters, combined with local geotechnical data. The following section explains the steps taken to determine model parameters. A summary of the used parameters is given in table S5.
As void ratio generally decreases with depth, we empirically derived a depth-dependent void ratio relationship based on a bulk analysis of almost 40,000 geotechnical samples from HCMC province (2) As sand is on average 20 times less compressible than clays and silts (Table 2.b in NNI, 2012), the for the sandy aquifers was subsequently corrected for this factor.
The recompression index ( ) was determined using local geotechnical data, following the linear relationship existing between the and indices (Gunduz and Arman, 2007). For the clayey aquitards, we found a / ratio of 5 based on the analysis of Mekong delta clays (Toan and Nu, 2013) and for the sandy aquifers a / ratio of 3, based on data from HCMC (Thoang and Giao, 2015).
Subsequently, the primary compression ratios, (compression ratio) and (recompression ratio), were derived from the and indices and the void ratio by: In principle, for a specific sediment, both compression ratios are independent of depth, however, the above approach introduces a decreasing trend with depth. Physically, this can be explained by a successive coarsening trend in lithology, which is the case for the Mekong delta (Fig. S4). The derived compression ratios for the shallow aquitards correspond to typical values associated with pure clays, and gradually decrease towards characteristic values for respectively weak sandy clay, sandy clay/loam and sandy loam for the deeper aquitards (Table 2.b in NNI, 2012).

Figure S4
Left: Depth-dependent compression ratio (CR) for aquitards in-line with a decreasing clay fraction with depth. Right: Lithoclass distribution of the 40% finest grained deposits in the Mekong delta, corresponding to the viscoplastic active part of the aquitards in the subsurface model. The coarsening trend in the lithoclasses with depth is clearly visible.
The coefficient of secondary compression ( ) for the aquitards was determined using a / ratio of 0.04 (Mesri and Godlewski, 1977;Ladd et al 1977). We assume no secondary compression taking place in the sandy aquifers, therefore is set to zero. For a similar reason, the sandy part in the aquitard layers (1/3 of the total layer thickness) was excluded from compression calculations.
As the abc model is based on natural (or Hencky) strain, the final step was to convert the , and to, respectively, the a, b and c. For low strain levels, which is the case in the model, the following relationships apply: / ln (10)

S5.2 Implication of the simplified groundwater-subsidence modelling approach
We used a one-way coupling approach to relate hydraulic head development and subsidence in the sense that we first model hydraulic head development with a conventional groundwater model, and then use the hydraulic heads as a function of space and time to drive the subsidence module. However, ideally or formally, a fully (two-way) coupled solution is required to ensure consistency between the strain (or compression) rates in the subsidence module on the one hand, and the storage rates in the groundwater model on the other. The simplified one-way approach was adopted because calibration of the fully-coupled model would demand excessive computation time.
To check the consistency between the independently calibrated specific storage values of the groundwater model and the compression behaviour of the subsidence module, we calculated elastic and 'virgin' specific storage values for the model aquitards from the compression (a or Cr or RR and b or Cc or CR) parameter values and estimates of the effective stress. The calculated elastic values were found to differ from the calibrated values by a factor ranging between 0.05 and 0.4, and the calculated 'virgin' values by a factor ranging between 0.3 and 2.2. The fact that these values are reasonably close to 1 indicates fair consistency. A model run in which the calibrated specific storage values of the aquitards were replaced by the compression-parameter based specific storage values (taking the mean value of the elastic and 'virgin' values) yielded slightly higher subsidence values (average value of cumulative delta-wide subsidence was ~2% higher). Thus although the parameters in the groundwater and subsidence model are not fully coupled, the potential impact on the modelled subsidence values falls well within the reported range of modelled subsidence associated with uncertainty in model parameters in general." Another uncertainty arises from the current aquitard discretization, which might not be refined enough in a vertical direction to correctly represent the delayed propagation of pressure decline in the less permeable layers. The majority of the model is behaving quasi-steady with very limited delays in pressure decline propagation. This implies that the current layer-discretization is sufficient for this delta-wide model, however it is not the case everywhere in the model and could still contribute to model uncertainty.
The percentages give the part of the total InSAR-measured subsidence reproduced by modelled groundwater extraction-induced subsidence. The calculation is illustrated by the following example for Fig. 9. Point A: 1.5/0.5 cm/yr, measured > modelled subsidence = 33% explained by the model; Point B: 1.5/1.5 cm/yr = measured equals modelled subsidence = 100% explained by the model; Point C: 0.5/1.5 cm/yr, measured < modelled subsidence = 100% explained by the model. The percentages in the figure are the average values of all cells.

Figure S7
Fit between modelled groundwater extraction-induced subsidence and InSAR-measured subsidence for the MKD (Erban et al 2014) upscaled to model resolution for the best estimate model. Linear fits (y=1x; modelled equals measured subsidence) is shown for the best estimate model (overconsolidation ratio (OCR) of 1.63) and, respectively, for the least and most conservative models with OCRs of 1.45 and 1.75 (data is not shown in this figure). The best estimate model has a rather low cross-correlation (r) of 0.28 which is to be expected as total InSAR measured subsidence and groundwater extraction induced subsidence do not correlate in absolute values (discussed in the main manuscript section: Subsidence in the Mekong Delta). A relative correlation between the two is present, shown by a Spearman's rank-order correlation (ρ) of 0.73. Subsidence rates are in average annual rates over the period 2006-2010.

S6.1 Resampling InSAR-measured subsidence rates
The spatial resolution of the InSAR-measured subsidence rates as presented by Erban et al (2014) (300 m 2 ) differs from the modelled subsidence cell size (1 km 2 ). To enable direct comparison, the InSAR-measured subsidence rates were resampled using a median bilinear interpolation creating a weighted value to match the model cell size and raster position.

S6.2 Calculation of spatial statistics on model results
When MKD is explicitly mentioned in the text, this means that model and InSAR measurements outside of the MKD, i.e. HCMC province and Cambodia, are not included in the analysis. As the extraction database does not provide a full spatial coverage of the delta, values of reported average subsidence caused by groundwater extraction are calculated only for parts of the MKD inside a 5 km radius of an extraction with extraction amounts exceeding 5 m 3 /km 2 /day ( Fig S4b). Average subsidence values reported for HCMC were calculated over the area for which >50 cm of total subsidence in the best estimate model was calculated (Fig. 6, main text).