A groundwater security model based on hydraulic turnover times and flow compartments.

Starting with a log-linear relationship between groundwater discharge per unit drainage area (Q/A b), hydraulic turnover time (t) and aquifer mobile storage (z), this study builds a groundwater security method at catchment scale. The method embeds previously published approaches to calculate Q/A b, t and z, and relies solely on stream flow discharges and watershed areas. The ability to build a method on a couple of variables is remarkable. The method recasts the calculated variables as aquifer security indicators (S Q, S t and S z), relating S Q with yield capacity, S t with self-depuration capacity and S z with resilience. Groundwater security is the weighted product of S Q, S t and S z. The method is validated with stream flow discharges and drainage areas concerning 294 hydrometric stations and their watersheds, located in continental Portugal. The results revealed a majority of moderately to highly secure watersheds, especially as regards S t (> 62%), while 7-10% were classified as very highly secured in general (S Q-S t-S z). The least secured basins are located in the more arid regions of continental Portugal (Northeast and south regions), as expected. The method can be easily transposed to any other region worldwide, with the necessary adaptions to regional climate, geological and topographic settings. • Compile stream flow discharge data and organize them as natural logarithms and logarithmic variations as function of time, to estimate Q, t and z; • Recast the Q, t and z values as S Q, S t and S z ratings, respectively, using the appropriate reclassification scales, and estimate watershed security levels, namely average security or customized (weighted) securities that highlight the contributions of Q/A b (watershed yield), t (aquifer's self-depuration capacity) or z (aquifer's resilience); • Use the results to draw illustrative diagrams and spatial distribution maps.

Specifications table

Method details

The method introduced in this study to assess groundwater security conceives the underground portion of river basins as compartments with different volumes and storage capacity, and relates groundwater discharge rates from the compartments with hydraulic turnover times. Within this general framework, the larger and more porous river basins that yield massive groundwater volumes every year and can sustain them for a long time in the absence of replenishment, are viewed as the most secure. This is because they can persistently respond to large water demands and are capable to neutralize potential contamination from the surface considering the long-lasting underground pathways that allow contaminant attenuation through filtration, adsorption or chemical and biological decay. On the opposite side are the small low-porosity watersheds that yield just small volumes of groundwater annually and can sustain them for just short periods. In this case, the reason to be the least secure is that these watersheds can barely respond to large water demands and are readily affected by infiltration of contaminated surface waters given the short turnovers that hamper self-depuration. In between these two end-member compartments, there will be an assortment of other catchments with intermediate dimensions and hydraulic properties, gradients and turnover times, and hence with intermediate groundwater security.

The method is composed of three modules. The first module was coined “groundwater compartment” module and established a general log-linear regression between groundwater discharge rate and hydraulic turnover time that allowed the definition of groundwater compartments as function of watershed volume and effective porosity (i.e., the catchment mobile storage). This relationship is the method's core, because it gathered the above-mentioned elements of groundwater security in a single equation. Modules 2 and 3 were developed because implementation of module 1 requires prior estimation of those elements. Thus, the second module, termed “turnover time”, estimated the hydraulic turnover time of a watershed from a previously published method of Pacheco [1], and improved that method to additionally assess the groundwater discharge rate. It also used the log-linear relationship of module 1 to estimate the mobile storage (z) from the calculated groundwater discharge rate (Q) and hydraulic turnover time (t). Finally, the third module recast Q, t and z as security indicators (S) of aquifer yield, groundwater quality and aquifer resilience, respectively, and determined an average security for the basin as product of SQ, St and Sz. This is why the module was called “security” module. The average security indicator can be customized to highlight the role of yield, quality or resilience, through power weighting of SQ, St or Sz, respectively. The three modules are described in detail in the next subsections.

Groundwater compartment module

The flow of groundwater (Q, L3T–1) through the section of an aquifer is proportional to flow velocity (v, LT–1) and effective section area (Ae = A × ne, L2, where ne is the aquifer's effective porosity):

On the other hand, velocity is a cinematic relationship between a path length (L, L) and the time required to travel along the path (t, T):

The replacement of Eqs. (2) in (1) gives:where Ve (L3) represents the effective aquifer's volume also known as mobile aquifer storage. If Q is the discharge of groundwater from a river basin that was assessed at time t, then Ve is the mobile catchment storage. If both terms of Eq. (3) are now divided by the basin's area (Ab, L2), the result is:where z = Ve/Ab (L) is the mobile catchment storage per unit basin area. If logarithms are finally taken from the left- and right-hand sides of Eq. (4), the result is:

A projection of (LT–1) as function of (T) describes a straight-line y = a + bx, with b = –1 and a = . Thus, it is possible to use this outcome to draw a general versus diagram composed of straight lines with slope –1 (b = –1) and intercept-y values pre-defined at certain values of . The line corresponding to z would represent a totally impermeable catchment (ne = 0) or the boundary between the surface and groundwater compartments (V = 0). Other lines can be given specific meanings. For example, a line drawn for z ≤ 1 m can be used to limit below the soil plus saprolite (altered rock layer) compartment whereas if drawn for 1 〈 z ≤ 10 m can be used to limit the shallow groundwater compartment. The deep groundwater compartment can be bounded below by the 10 < z ≤ 100 m lines and finally the very deep groundwater compartment by the z 〉 100 m line. These boundaries are not universal but are acceptable for Portuguese watersheds considering results from various studies [2], [3], [4], [5]. Other geologic, topographic and climate settings may change these limits, but users can adapt the boundaries to their own settings. The diagram can be further divided by vertical lines representing pre-defined hydraulic turnover times. Thus, within the z compartments, there will be regions of fast (e.g., t < 1 year) to slow (e.g., t > 100 year) flow, representing watersheds of similar storage capacity but with different hydraulic conductivity and/or gradient. Taken altogether, the limits imposed to z and t, coupled with the range of Q/Ab values, define sectors in the diagram that can be interpreted from the standpoint of groundwater security. In general, the larger the Q/Ab, z and t the best for security. Larger values of Q/Ab ensure improved capacity to attain water demand for activities while larger values of z warrant resilience against hydrologic drought. Moreover, larger values of t safeguard longer contact between dissolved pollutants and aquifer materials that allow more efficient contaminant attenuation through filtration, adsorption or decay processes. Thus, starting from a simple log-linear regression (Eq. (5)), the proposed method is capable to assess water security from quality as well as quantity viewpoints.

Having defined the groundwater flow compartments, as function of Q/Ab, z and t, the diagram can be populated with and values relative to actual watersheds. The plot of a real watershed over the compartments allows to envisage which type the watershed represents. However, before doing that, a method is required to estimate the groundwater discharge and the hydraulic turnover time at catchment scale. The present study resorted to a model developed by Pacheco and published in 2015 [1] that estimates t as function of stream flow discharge, which will now be improved to estimate groundwater discharge as well. A detailed description of the method is presented in the next subsection.

Hydraulic turnover time module

Similarly to Eq. (4), the hydraulic turnover time of a watershed has been equated to [6]:where t (T) is the time, Vne (L3) is the mobile catchment storage, with V (L3) equated to the watershed volume and ne (dimensionless) to its concomitant effective porosity, and Q (L3/T–1) is the average groundwater discharge. A direct use of Eq. (6) to estimate the hydraulic turnover time implies the prior assessment of its components, but Pacheco [1] reduced the required data to Q when the equation is combined with a formula for ne deduced from a recession flow method, namely the Brutsaert method [7], [8], [9].

The Brutsaert method used the so-called Boussinesq equation established in 1903 to describe the drainage from an ideal unconfined rectangular aquifer bounded below by a horizontal impermeable layer and flowing laterally into a water channel. The solution for that equation has the general form of a power function:where Q (L3/T–1) represents the groundwater discharge and t (T) is the time, while a and b represent hydraulic coefficients. The a coefficient relates with the groundwater reservoir's characteristics and b with the stream flow regime, namely the short-time or high-flow (b = 3) and the long-time or low-flow (b = 1) regimes. In the work published in 1998, Brutsaert and Lopez [9] derived the following solution for the short-time regime:where K (LT–1) is the hydraulic conductivity, ne (dimensionless) the effective porosity, z (L) the aquifer thickness and l (L) the length of upstream channels intercepting groundwater flow. The long-time regime is adequately described by the so-called linear solution of Boussinesq published in 1903:where Ab (L2) is the upland drainage area. In his work of 2015, Pacheco noted that when Eqs. (8a) and (8b) are combined and z is equated to V/Ab, the effective porosity becomes written as:where ai represents the value of a when b = i (1 or 3). Besides, He further noted that by combining Eq. (9) with Eq. (1), the estimation of hydraulic turnover time simplifies to:meaning a formulation solely dependent on average groundwater discharge (Q) and flow regimes (a1 and a3). According to the Brutsaert method, the values of a1 and a3 can be read in a scatter plot ln(DQt/Dt) versus ln(Qt), where Qt is the stream flow discharge measured in a hydrometric station located at the outlet of a catchment. In that plot, the lower envelope to the scatter points is represented by two straight lines, one with a slope b = 1, and the other with a slope b = 3, and the y-values where these lines intercept ln(Qt) = 0 are the parameters a1 and a3.

In the present study, the average groundwater discharge (Q) is also determined from the scatter plot ln(DQt/Dt) versus ln(Qt). In that diagram, the lines of slope 1 and 3 intersect each other at Qmax, which represents the maximum groundwater discharge in the low-flow regime because the intersection point separates the low-flow regime (only groundwater discharge) from the high-flow regime (groundwater + surface water discharge). In the present study, Qmax is used as proxy to the maximum possible groundwater discharge. To estimate Qmax, the user draws a vertical line through the intersection point of slope 1 and 3 lines and reads ln(Qmax) at the independent variable axis (X-axis). On the other hand, the scatter point plotted most to the left in the low-flow regime region represents the lowest possible groundwater discharge in the available stream flow record (Qmin). In this case, drawing a vertical line through this point allows the reader to estimate ln(Qmin) at the intersection of this line with the independent variable axis. Finally, the (geometric) average groundwater discharge will be given by:

After the determination of Q and t, z can be estimated using Eq. (5).

Security module

As mentioned above, variables Q/Ab, t and z can be recast as indicators of groundwater security. Variable Q/Ab describes security from a quantity (yield) standpoint, as larger yields ensure a better response to water demand and especially peaks of demand. Variable t describes security from a quality viewpoint because of its relationship with pollution attenuation. And finally, variable z looks at security from the side of aquifer resilience, because a larger mobile storage allows a better adaption of water supply systems to variations of aquifer replenishment through infiltration of precipitation, namely to long periods of meteorologic and hydrologic drought.

The security module recasts Q/Ab, t and z as security indicators (SQ, St and Sz, respectively) through reclassification. Considering the log-linear relationship between the three variables (Eq. (5)), reclassification will be also log-linear. The security indicators increase as function of increasing values of their parent variables, because yield, quality and resilience security increase as Q/Ab, t and z increase. The security classes cannot be defined universally because they are likely dependent on regional settings such as climate, geology or topography. For the present study, adapted to the pilot application in continental Portugal, the security classes were defined as depicted in Table 1.

Table 1

Security classes of Q/Ab, t and z. They were adapted to the pilot application that spans a large number of Portuguese watersheds.

Security class	Q/Ab range (m/year)	t range (year)	z range (m)
1	< 0.01	< 10	< 1
2	0.01–0.1	10–100	1–10
3	0.1–0.5	100–500	10–50
4	0.5–1	500–1000	50–100
5	> 1	> 1000	> 100

The average security of a river basin is defined by a weighted product of partial securities (SQ, St and Sz), i.e.:where wq, wt and wz are the weights (importance) attributed to the SQ, St and Sz indicators, respectively. By default, the weights are all equated to 1 (all indicators are equally important), but they can be customized to raise the importance of a specific indicator. In that case, the weight of an indicator can be set to values between 1 and 3, providing that the other weights are set to values between 0 and 1, and that the sum of weights is always 3, i.e.

Setting the sum of weights fixed to a constant value keeps the range of S values relatively uniform and hence allows a direct comparison among weighted and unweighted security diagrams or maps. If the weight of an indicator is set to 3 (and the other weights to 0), then the resulting diagram or map will represent groundwater security from a single standpoint (yield, quality or resilience).

The calculated S values (Eq. (12)), if weights are all 1, range from 1 to 125. The final step of security module is to set up a qualitative scale that levels groundwater security between very low and exceptional, hinged on five classes of S. The criterion used to define the class boundaries was: (1) class 1 is defined when all indicators are 1 (i.e., Sq = 1 and St = 1 and Sz = 1); (2) the class is set to i + 1 if the calculated S is larger than a reference case where two indicators are equal to i and one indicator is equal to i + 1. Based on this criterion, the levels and ranges of groundwater security are depicted in Table 2.

Table 2

Levels of groundwater security and corresponding ranges of S as determined by Eq. (12).

Level of security	S range
Very low	1
Low	1–2
Moderate	2–12
High	12–36
Very high	36–80
Exceptional	80–125

Data and software requirements

The groundwater security method is minimalist as regards data requirements. It fully operates with a dataset of daily stream flow discharges measured at the outlet of target watersheds, complemented with watershed boundaries and areas. In general, the records of stream flow discharge can be downloaded from the websites of public water resource management institutions, and saved as spreadsheets. Sometimes, these institutions also provide shapefiles / geodatabases with the delineation and geometric characterization of watersheds located upstream the sites where the stream flow discharges were measured. If not available, the stream flow discharges can be measured on site with appropriate hydrometric stations, and the watershed boundaries interpreted from Digital Elevation models using conventional terrain modeling tools embedded in GIS software.

The three modules comprising the groundwater security method can all be implemented in Microsoft Excel software using solely the stream flow discharge measurements and watershed areas (e.g., Fig. 1a,b below). The results (e.g., security indicators) can then be used in more sophisticated statistical computer packages for appealing representations. Here, the STATISTICA software of Statsoft (https://www.statistica.com/en/) was used to draw contour plots (e.g., Fig. 2a–d) and histograms (Fig. 3a–d), whereas the ArcGIS Pro-of ESRI (https://www.esri-portugal.pt/pt-pt/arcgis/produtos/arcgis-pro/overview) was used to draw spatial distribution maps (Fig. 4a–d).

Fig. 1

a – Example of how to apply the Brutsaert method to stream flow discharge data. The example refers to station 03 J/02H. The raw data and full implementation of the method can be consulted in the Supplementary Materials (Spreadsheet 2). 1b – Distribution of 294 Portuguese watersheds (blue circles) as function of groundwater discharge (Q/Ab), hydraulic turnover time (t) and watershed mobile storage (z).

Fig. 2

Distribution of groundwater security (shaded areas) as function of groundwater discharge (Q/Ab), hydraulic turnover time (t) and catchment mobile storage (z): (a) Saverage – unweighted security; (b) Squality – groundwater security highlighting the role of hydraulic turnover time; (c) Syield - groundwater security highlighting the role of Q/Ab; (d) Sresilience - groundwater security highlighting the role of z.

Fig. 3

Groundwater security histograms: (a) Saverage; (b) Squality; (c) Syield; (d) Sresilience. Additional information in the caption of Fig. 2.

Fig. 4

Groundwater security of continental Portugal watersheds, as estimated with the proposed method: average(a) and highlighting groundwater quality protection (b), aquifer yield (c) and aquifer resilience (d).

Strong points and gaps

Strong point 1: novelty. The method has the capacity of unraveling three important viewpoints of groundwater security using a single equation (Eq. (5)), namely aquifer yield capacity, aquifer self-depuration capacity and aquifer mobile storage (resilience to prolonged drought periods), which is remarkable. Strong point 2: minimal data requirements. The method relies on a couple of readily available variables: daily stream flow discharge and watershed area. Gap: parameter estimation dependent on user experience. The subjectivity of determining a1, a3, Qmin and Qmax from Fig. 1a is a source of uncertainty.

Method validation

Data preparation and method implementation

The method is exemplified with streamflow discharges and drainage areas relative to 294 watersheds located in continental Portugal, spanning an amble range (several orders of magnitude) of these values: 0.02 < stream flow discharge (m3/s) < 517.13; and 2.6 × 106 < drainage area (m2) < 9.7 × 1010. The data were retrieved from the Portuguese System for Information on Water Resources, at the URL http://snirh.apambiente.pt. The Supplementary Materials (Spreadsheet 1) contain the list of basins and of hydrometric stations located at the basin outlets, namely information on basin's name and drainage area, as well as on station's code, name, geographic location (latitude, longitude, altitude) and average stream flow discharge. Spreadsheet 2 can be consulted to see how the Brustsaert method was applied to one of the hydrometric stations (03 J/02H), namely: (1) how the values of a1 and a3 were estimated from the ln(DQt/Dt) versus ln(Qt) diagram; (2) how the values of ln(Qmax) and ln(Qmin) were estimated from the same graph; and (3) how the previous values were coupled with the average Q value (Eq. (11)) to estimate t using Eq. (10). The stream flow discharges depicted in Spreadsheet 2 refer to daily records, which were summarized as monthly averages before being used with the Brutsaert method. The scatter diagram plotting ln(DQt/Dt) versus ln(Qt) is also represented in Spreadsheet 2 and reproduced in Fig. 1a for illustrative purposes. The DQt/Dt quotient was estimated as (Qti+1 − Qti)/(t+1 − t), where subscript i refers to a month. The graphic's abscissa was estimated as ln[(Q+1 + Q)/2]. The dashed lines were drawn with slope 1 and slope 3 just below the cloud of points. It is worth noting the haziness of drawing these lines, because it is frequent that some points scatter significantly below the main cloud. In those cases, the dashed lines can be drawn so 90% of the points are located above them [10,11]. A vertical orange line was passed through the intersection of slope 1 and slope 3 lines to obtain ln(Qmax). Another vertical line was passed through the lowest continuous set of ln(Qt) values to obtain ln(Qmin). Two points at the left of this line were not considered in this assessment because they were assumed outliers. The intercept-y of the dashed lines in Fig. 1a are ln(a1) = −20 and ln(a3) = −25 which, for an average groundwater discharge of Q = 2.34 m3/s (or 7.38 × 107 m3/year), implies a hydraulic turnover time of t = 158.6 years for station 03 J/02H (River Vez basin). By drawing plots similar to Fig. 1a, but relative to the other 293 hydrometric stations listed in the Supplementary Materials, hydraulic turnover times and average groundwater discharges were calculated for all studied basins. The results are listed in Spreadsheet 3 of Supplementary Materials and range from: t = 1.39 to 3185.79 years; Q/Ab = 0.003 to 11.32 m/year. It is worth recalling here that turnover times determined by Pacheco in 2015 [1] were based on average Qt and not on average Q and hence are likely underestimated.

Having estimated the hydraulic turnover times for all river basins, the 294 Q/Ab versus t points were plotted in the groundwater security diagram (Fig. 1b). The studied basins span various groundwater compartments (various orders of z), which is not surprising given the range of drainage areas that span > 4 orders of magnitude (2.6 × 106 < drainage area (m2) < 9.7 × 1010).

A brief summary of results in the tested area

The groundwater security of watersheds is represented in Fig. 2a–d. The allocation of security ratings to the watersheds is in keeping with the classes of Table 1 (security indicators) and Table 2 (security levels as determined using Eq. (12)). The allocation of SQ, St and Sz, as well as the calculation of S, are provided as Supplementary Materials (Spreadsheet 4). Fig. 2a describes unweighted security while Fig. 2b–d result from a customization of weights. In the case of Fig. 2b, wt = 1.5 and wq = wz = 0.75, meaning that this figure highlights the contribution of hydraulic turnover time (protection of groundwater quality) to groundwater security. Fig. 2c and d were drawn with the purpose to highlight the roles of yield and resilience, respectively, and, in those cases, the weights were set up as follows: wq = 1.5 and wt = wz = 0.75 (Fig. 2c); wz = 1.5 and wq = wt = 0.75 (Fig. 2d). In all these figures, the vast majority of watersheds plot where security scores range from 2 to 36, meaning that they are moderately to highly secured. There are even a significant number of watersheds with very high security (points plotted over the red-shaded areas). The histograms of security, with identification of counts and percentages of count, are presented in Fig. 3a–d, as complement to the diagrammatic representation of Fig. 2a–d. These figures point to 7–10% of very high security watersheds and > 60% watersheds highly secured for groundwater quality.

The spatial distribution of groundwater security is displayed in Fig. 4 for the average (a), quality (b), yield (c) and resilience (d) scenarios. It is evident the concentration of moderately secure basins in the more arid regions of continental Portugal (Northeast and South, where precipitation is low), and the general highly secure basins concerning quality.

Supplementary material The Supplementary Materials are composed of four spreadsheets: (1) list of watersheds used to validate the groundwater security method. The spreadsheet contains information on basin's name and drainage area; the code, name and location information about the hydrometric station located at the basin's outlet; (2) exemplification of how to apply the Brutsaert method to one hydrometric station (Fig. 1a); (3) summary of hydraulic turnover times and projection of Q/Ab versus t values on the groundwater security diagram (Fig. 1b); (4) Base data and calculation of security indicators (SQ, St and Sz) and security levels (Saverage, Squality, Syield and Sresilience) used to produce Fig. 2 (diagrammatic representation of security levels as function of groundwater discharge, hydraulic turnover time and catchment mobile storage), Fig. 3 (histograms of security levels) and Fig. 4 (spatial distribution of security levels).

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Subject Area;	Environmental Science
More specific subject area;	Groundwater hydrology
Method name;	A groundwater security model based on hydraulic turnover times and flow compartments
Name and reference of original method;	Pacheco, F.A.L., 2015. Regional groundwater flow in hard rocks. Science of the Total Environment 506–507, 182–195.
Resource availability;	Besides the Supplementary Materials provided with this submission, the authors are committed to provide additional data at the request of an interested scientist