Quantifying Groundwater Flow Using Water Budgets and Multiple Conservative Tracers

 S.V. Smith, V. Camacho-Ibar, J. Herrera-Silveira, D. Valdes,

L. David, M. Merino and R.W. Buddemeier

 Introduction

 One of the most difficult estimates for C, N and P biogeochemical budgeting in many systems is the input of groundwater (V_G) and its dissolved constituents into coastal ecosystems.  In regions such as the Yucatán Peninsula of south-eastern México, where groundwater flows are known to be significant and surface runoff is near 0, this determination becomes critical (see, for example, Perry and Velasquez-Oliman 1996; Back and Hanshaw 1970; Hanshaw and Back 1980).  In this particular environment, Hanshaw and Back estimated that the 1,100 km northern portion of the peninsula has an average groundwater outflow to the ocean of approximately 8.6 x 10⁶ m³ km^-1 yr^-1.  Those authors report some evidence for spatially heterogeneous distribution of this outflow, as would be anticipated from the very heterogeneous distribution of sinkholes, or cenotes (e.g., Perry and Velasquez-Oliman 1996).

 A geochemical approach to the estimation of groundwater flow is the use of alternative geochemical tracers.  Examples include radium (Moore 1996a and b), radon (Burnett et al. 1996), and methane (Chanton et al. 1996).  Data provided by several participants to this workshop provided hints that silicate might be used as a quasi-conservative tracer of groundwater inputs to some coastal lagoons around the Yucatán Peninsula.  D.R. Corbett (personal communication) has noted that groundwater silicate appears to be nearly conservative in the carbonate terrain of Florida Bay.  In the case of Yucatán, silicate levels are elevated and not highly variable in groundwater and surface waters across much of the peninsula (Herrera-Silveira et al. 1998; Herrera Silveira 1999), and several of the coastal lagoons discussed in this report show elevated Si concentrations (e.g., Celestún; Herrera et al. this report).   Celestún Lagoon provides particular insight.   Mixing diagrams suggest that the silicate distribution in this system with known groundwater discharge may be distributed approximately conservatively with respect to salinity (Herrera-Silveira 1995; Herrera-Silveira and Ramírez-Ramírez 1998).   Because the atmospheric term (precipitation minus evaporation) in many of the Yucatán lagoons is significant relative to groundwater flow, mixing diagrams of silicate versus salinity will not in general be linear; there are three end-member water masses, not two (see Boyle et al. 1974). 

Chelem Lagoon and Ria Lagartos are of interest in this context, because these systems have elevated silicate (hence, the suggestion of groundwater input) even though they are hypersaline systems (Valdés and Real 1998; Valdés this report).  There is, moreover, precedent for using Si as a hydrological tracer (Kennedy et al. 1986; Wels et al. 1991).  

 The Nichupté Lagoonal System also provides insight.  Systems such as this one appear to receive significant freshwater inflow from local runoff which occurs during rainfall events (Merino et al. 1990), and this complication must be considered.

Theoretical analysis

With this background in mind, we have expanded the calculations laid out in Gordon et al.  (1996), to the case where groundwater might have a unique signature which would distinguish it from rainwater.  Conservation of water in the system is described as follows:

                                                                                            (1)

 where  dV_sys /dt is the change of volume in the system over time; the V's denote volume fluxes, with the subscripts Q, P, E, G, and O representing river discharge, precipitation, evaporation, groundwater, and other freshwater sources (e.g., sewage), respectively.  The subscript R is the ‘residual flow’ necessary to conserve mass and is treated as the unknown.  It is convenient, for the discussion to follow, to consider that V_PE represents the net of precipitation minus evaporation.  Moreover, it is convenient (although by no means necessary) to treat the system as at steady state, that is dV_sys/dt = 0.  With these simplifications and assumptions, the equation can be solved for V_R:

                                                                                            (1a)

 If salt is conserved, then a similar equation can be written for the conservation of salt (S).  The equation will have an additional term (V_X), to describe the mixing of water between the system of interest and the ocean (ocn):

                   (2)

 In this equation, S_sys and S_ocn are the salinity values for the system and oceanic boxes; S_R, the ‘residual salinity’, is taken to be the salinity at the boundary between the system and the ocean (i.e., the average of S_ocn and S_sys).  Some of the terms (e.g., salinity of the precipitation - evaporation and the term for other flow) can usually be treated as near 0 and are dropped out of the analysis here in order to simplify the equations; they can be re-inserted in systems where they might be quantitatively significant.  One clear test of ‘significance’ is to evaluate both these simple versions of the equations and the more complete equations - even if with hypothetical data.  A change of less than 25% in the estimate of unknowns probably lies within the range of uncertainty in the known quantities.  We can consider this equation at steady state (that is, VdS/dt = 0), combine it with (1a), and solve for V_X:

                                                      (3)

 If all the terms on the right side of the equation are known, then obviously an estimate for V_X can be derived; this is the ‘standard LOICZ procedure’.  Consider the case where one term (V_G) is not known, but for which a second conservative tracer (let us assume that tracer is silicate, Si) is known.   An equation exactly analogous to (3) can be written:

                                                           (4)

 Equations (3) and (4) can be combined to give an estimate of V_G:

                      (5)

Results and Discussion

There are several caveats to the use of equation (5).  The most obvious ones are that the salinity and silicate content of the groundwater are known, and that the silicate concentration is conservative with respect to salinity.   Slight nonconservative behavior of the silicate would ordinarily be expected to introduce a relatively small error in the calculations.   Because these systems receive little or no surface-water discharge, which might have high concentrations of diatoms (with relatively soluble SiO₂ frustules), reactive forms of particulate silicate are not being supplied to these systems.  Indeed, the classical paper by Boyle et al. (1974) reviewing the use of mixing diagrams to assess the chemical mass-balance of estuaries concluded (p. 1724): “...in no case has it been proved unambiguously that silica exhibits non-conservative behavior in estuarine mixing.”  Moreover, because groundwater and lagoonal silicate concentrations are very high compared to concentrations of dissolved inorganic N and P, large deviations from conservative behavior are not expected in these systems. 

A second set of considerations may actually be more important.  Not surprisingly, V_G is calculated as a volume flux scaled to the other freshwater input terms - V_Q, and V_P (i.e., the net of V_P and V_E) as the equation has been simplified and formulated.  The value for V_G is therefore only as good as the estimates of these other flux terms.  Moreover, the calculation of V_G as formulated is actually a calculation of the flux of water with high silicate (and usually low salinity).  V_Q, as well as V_G, is likely to fit that profile.  Therefore, if V_Q is either poorly constrained or large relative to V_G, then the calculation will not be robust.  Of course some other tracer might differ between river-water and groundwater and could be substituted for silicate.  In the case of the northern Yucatán Peninsula, with virtually no river flow and high groundwater flow, the equation generally appears robust.

A third caveat, which we have learned by application of this equation to examples in this report, is that the estimate of V_G is quite variably sensitive to the estimated values of salinity and silicate in the groundwater.  In some instances, the calculation is sensitive to one of these variables, sometimes to the other, sometimes to both, and sometimes to neither.  The sensitivity is largely reflecting regions where the denominator of equation (5) is close to 0 and responsive to slight variations in these two variables.  We have found it convenient to create a spreadsheet matrix with a range of salinity and silicate values (for most cases in the Yucatán systems, salinity ranging from 0 to 12 psu in steps of 1 psu, and silicate ranging from 100 to 500 mM, in steps of 50 mM, appears appropriate).  In effect, for the observed characteristics of salinity and silicate concentrations in the lagoon, this matrix is a ‘sensitivity map’ in salinity-silicate space.   We can then look in the matrix to see where the estimated value falls with respect to sensitivity to these two variables.  In cases with the denominator near 0, slight variations can make the estimated value for V_G become either very large or negative or both.   

An abbreviated version of this matrix is shown as Table 1 and graphically as Figure 1, for Celestún Lagoon.  In this example, estimated groundwater flux is in a region of the matrix that is moderately sensitive to uncertainty in salinity higher than the estimated value and very sensitive to higher silicate concentrations.  Lower values of either salinity or silicate do not dramatically alter the calculations. 

Versions of equation (5) are used in several of the case studies given in this report: Celestún Lagoon and Dzilam Lagoon (both budgets by Herrera-Silveira et al.), and Chelem Lagoon and Ria Lagartos (Valdez) (Table 2).  All of these systems are located in the state of Yucatán, along the north and north-west portion of the Yucatán Peninsula.  None of these systems has significant river inflow; all have evaporation in excess of precipitation; Celestún and Dzilam both have salinity below oceanic, even though evaporation exceeds rainfall in the region.  Chelem and Ria Lagartos are both hypersaline throughout most of their extent, although salinity at the mouths of these systems is slightly lower than coastal seawater (Table 3).   The presence of depressed salinity at the mouths of these systems, even though they are in net evaporative regions with no significant river flow, is proof that a low-salinity source (groundwater) must be important in the water budgets.  All four of these systems show elevated silicate levels in the lagoon waters.  The range of estimated V_G for these systems taken as whole units is about 1 to 4x10⁶ m³ km^-1 yr^-1 (Table 4), with portions of the systems showing locally much higher rates (see individual nutrient budgets in main body of this report). 

A system for which the calculations did not initially seem to work was the Nichupté Lagoonal System, Quintana Roo (Merino).  There, the initially calculated groundwater fluxes were negative (Table 5A).   We recognise that negative groundwater flux (i.e., saline intrusion into the aquifer) does occur in some locations.  Indeed, that is a significant problem in many areas of México where groundwater exploitation exceeds recharge.  This is not the case in most of Yucatán, because of the large volume of recharge and relatively low utilisation rates.  Moreover, there are known springs in Nichupté. 

After examination of the salinity-silicate sensitivity matrix, we think that the problem lies with local surface flow which is not adequately accounted for in the water budget.  For this system, V_P - V_E for the analyzed period was -21x10⁶ m³ yr^-1.  A positive flux for non-groundwater freshwater inflow would reverse the sign of the estimated groundwater flux.  This interpretation is consistent with the analysis by Merino et al. (1990).  Those authors observed that the wetlands immediately adjacent to Nichupté cover an area approximately equal in size to the lagoon, and that runoff from a significant fraction of this wetland area is apparently important.  Those authors estimated that local runoff during rainfall events might deliver between two-thirds and all of that rainfall directly to Nichupté.  It can be assumed that this local runoff would be low in salinity and probably would have had inadequate time to have elevated silicate concentrations.  For the period in question, adding 67% of the rainfall as local runoff would be equivalent to adding 44x10⁶ m³ yr^-1 of additional fresh water.  When this local runoff is added (Table 5B), the resultant estimates of groundwater flux become positive.

We also tried the use of the equation for Terminos Lagoon, Campeche (David), on the southwest portion of the Yucatán Peninsula, without success.  Calculated groundwater flow was clearly far too large to be physically reasonable, although it was still a small quantity in comparison to river flow.  In that instance, the dominance of freshwater inflow by rivers (also high in silicate and low in salinity) precludes the ready use of this equation.  

Conclusions

We conclude that the use of silicate as a second ‘conservative tracer’ seems to work as an estimator of groundwater flow for much of the northern Yucatan Peninsula.  Undoubtedly this technique might be further adjusted, especially with site-specific data on groundwater composition.  Moreover, it is clear that specific considerations such as local runoff should be taken into account in the water budget.  Finally, domination of the water budget by river flow, which is likely to have a silicate concentration similar to that of groundwater, will compromise this approach.

The water fluxes associated with groundwater in the northern Yucatan Peninsula are significant to both the water and nutrient (especially nitrogen) budgets of the lagoons (main body of report).  The work by Corbett et al. (in press) in Florida Bay underscores the potential importance of groundwater in the nutrient budgets of such carbonate terraines with high groundwater flow and low surface flow.

Although the mean flow rates at the scale of the individual systems in the Yucatan Peninsula appear to be well below the regional estimate of Hanshaw and Back (1980), we believe that the general pattern is consistent with their analysis; it seems likely that much of the Peninsula does, indeed, have low groundwater flow rates, and that small regions account for a significant proportion of the total flow for the entire region. 

Table 1.   Salinity--silicate sensitivity matrix for Celestún Lagoon.   Rainfall minus evaporation for this system is -17x10⁶ m³ yr^-1; lagoon mouth and oceanic salinity and silicate values are given in Table 3.  As summarised in Table 4, the estimated groundwater flow (at groundwater salinity and silicate concentrations of 9 psu and 244 mM, respectively) is 51x10⁶ m³ yr^-1.  Figure 1 illustrates this same matrix graphically.

Table 2.  Physical dimensions and estimated rainfall minus evaporation data for four Yucatán coastal lagoons.  

Table 3.  Estimated water composition for groundwater, water at the mouth, and open coastal seawater, for the four lagoons listed in Table 2.  In the case of Celestún, annual average data are reported here; the text in the main body of the report uses seasonal data.  For Chelem and Lagartos, groundwater salinity and silicate data are estimated from Herrera-Silveira et al.  (1998). 

Table 4.  Estimated groundwater fluxes for the four lagoon systems listed in Table 2, based on data in Tables 2 and 3, and solution of Equation 5.  Rounding differences and seasonal versus annual data result in slight discrepancies between the data reported here and that in the main body of the report. 

Table 5.   Salinity--silicate sensitivity matrices for Nichupté Lagoonal System.   Part A is calculated with V_P-V_E = -17x10⁶ m³ yr^-1 and without local runoff.  Lagoon salinity and silicate values are 27.7 psu and 7 mM, respectively; oceanic values are 31.7 and 2.  Note that over this apparently reasonable range of groundwater salinity and silicate values, the estimated groundwater flux is consistently negative.  Part B repeats the calculation but adds 44 x 10⁶ m³ yr^-1 of local runoff, as adapted from Merino et al. (1990).