Rebuttal to “The case of the Biscayne Bay and aquifer near Miami, Florida: density-driven flow of seawater or gravitationally driven discharge of deep saline groundwater?” by Weyer (Environ Earth Sci 2018, 77:1–16)

A recent paper by Weyer (Environ Earth Sci 2018, 77:1–16) challenges the widely accepted interpretation of groundwater heads and salinities in the coastal Biscayne aquifer near Miami, Florida, USA. Weyer (2018) suggests that the body of saltwater that underlies fresh groundwater just inland of the coast is not a recirculating wedge of seawater, but results instead from upward migration of deep saline groundwater driven by regional flow. Perhaps more significantly, Weyer (2018) also asserts that established hydrologic theory is fundamentally incorrect with respect to buoyancy. Instead of acting along the direction of gravity (that is, vertically), Weyer (2018) claims, buoyancy acts instead along the direction of the pressure gradient. As a result, Weyer (2018) considers currently available density-dependent groundwater flow and transport modeling codes, and the analyses based on them, to be in error. In this rebuttal, we clarify the inaccuracies in the main points of Weyer’s (2018) paper. First, we explain that Weyer (2018) has misinterpreted observed equivalent freshwater heads in the Biscayne aquifer and that his alternative hypothesis concerning the source of the saltwater does not explain the observed salinities. Then, we review the established theory of buoyancy to identify the problem with Weyer’s (2018) alternative theory. Finally, we present theory and cite successful benchmark simulations to affirm the suitability of currently available codes for modeling density-dependent groundwater flow and transport.


Introduction
Lateral intrusion of seawater is a common occurrence in coastal aquifers around the globe (for example, Barlow and Reichard 2010;Bocanegra et al. 2010;Custodio 2010). In a study of saltwater intrusion into the Biscayne aquifer near Miami, Florida, USA, Cooper et al. (1964) were among the first to present detailed field observations and interpretation of circulatory flow within a "saltwater wedge," in which seawater moves inland at depth, reverses direction, and mixes with overlying freshwater before discharging at the coast. Their widely accepted interpretation of flow within the saltwater wedge is based on field measurements of salinity and saltwater head, calculations of equivalent freshwater head, and a dye tracer test. Numerous numerical studies of saltwater wedges (for example , Smith 2004;Abarca et al. 2007a, b) have since examined the effects of discharge of freshwater to the sea, the density difference between freshwater and saltwater (that is, buoyancy), hydrodynamic dispersion at the freshwater-saltwater interface, and aquifer geometry on the type of circulatory flow identified by Cooper et al. (1964).
In a recent paper published in Environmental Earth Sciences titled "The case of the Biscayne Bay and aquifer near Miami, Florida: density-driven flow of seawater or gravitationally driven discharge of deep saline groundwater?" Weyer (2018) contends that "[t]he salt wedge interpretation forwarded by Kohout (1964, in Cooper et al. 1964) constitutes a severe misunderstanding and misinterpretation of field data" because it was formulated "without taking nonchemical field data into consideration," despite the head and tracer observations presented and discussed by Kohout (1964) in the same publication. In fact, Weyer (2018) uses equivalent freshwater heads plotted by Kohout (1964) to reinterpret saltwater flow directions and concludes that "the so-called seawater wedge was actually caused by discharging deep saline groundwater driven by regional gravitational groundwater flow systems," rather than by density-driven flow of seawater.
Even more significant are the implications of Weyer's (2018) alternative conceptualization of density-dependent groundwater flow, which draws on Hubbert's (1940Hubbert's ( , 1953 landmark work on potentials and driving forces. Based on Hubbert's (1953) analysis of the effect of density on the force per unit mass acting on a fluid body, Weyer (2018) concludes that rather than act vertically (along the direction of gravity), buoyancy induces a driving force "along the direction of the pressure potential forces," that is, along the direction of the pressure gradient. As this view of buoyancy runs counter to established fluid dynamical theory, Weyer (2018) asserts that "[h]ere a paradigm shift occurs which contradicts the principles of 'density-driven flow' in any system with hydrodynamic boundary conditions" and that "the new understanding needs to be included in numerical modeling and proper practical water supply management world-wide." In this rebuttal, we address the inaccuracies in Weyer's (2018) interpretation of saltwater intrusion into the Biscayne aquifer near Miami and his view of density-dependent flow in general. First, we discuss Weyer's (2018) misinterpretation of observed equivalent freshwater heads in the Biscayne aquifer and critique his hypothesis concerning the source of the saltwater. Then, we review Hubbert's (1953) analysis of pressure and gravitational forces to explain the problem with Weyer's (2018) conceptualization of buoyancy. Finally, we verify that the forms of Darcy's Law employed in commonly used density-dependent groundwater flow and transport simulation codes are consistent with the classical theory of Hubbert (1940).

Misinterpretation of equivalent freshwater heads in the Biscayne aquifer
Density-dependent groundwater flow problems are often cast in terms of equivalent freshwater head (for example, Guo and Langevin 2002). If a piezometer is open to a given point within an aquifer, the equivalent freshwater head is the height of a column of freshwater in the piezometer that would be needed to balance the pressure at that point in the aquifer. In an aquifer filled with hydrostatic freshwater, the equivalent freshwater head is uniform with depth (and equal to the hydraulic head). In an aquifer filled with hydrostatic seawater, pressure increases more rapidly with depth than in the case of the less-dense freshwater. Thus, the equivalent freshwater head in a hydrostatic column of seawater increases with depth, and yet, no vertical flow occurs. This illustrates the well-known fact that a vertical gradient in equivalent freshwater head alone is not necessarily indicative of vertical flow. Kohout (1964) correctly interprets the contours of equivalent freshwater head, which he calls equipotential lines, within both the freshwater and saltwater parts of the Biscayne aquifer adjacent to Biscayne Bay. Within the upper, freshwater part of the aquifer, "flow lines must be nearly perpendicular to these equipotential lines," which implies "a seaward movement of freshwater." Within the underlying, saltwater part of the aquifer, equivalent freshwater head alone is not sufficient to fully determine the flow direction. However, under the assumption of a homogeneous aquifer, the horizontal component of the equivalent freshwater head gradient does indicate the horizontal direction of saltwater flow. Kohout (1964) correctly bases his conclusions regarding the direction of saltwater flow on the direction in which the equivalent freshwater head contours are inclined, which indicates the direction of the horizontal gradient.
In contrast, Weyer (2018) incorrectly equates the twodimensional direction of flow (which includes both the horizontal and vertical components) with the direction of the equivalent freshwater head gradient throughout the aquifer. This leads to unorthodox conclusions regarding the flow field and the source of the saltwater. Although Weyer (2018) concedes that the flow directions thus obtained are "approximate," they are inaccurate enough so as to be misleading. Weyer (2018) contends that his conceptual model of regional groundwater flow is supported by the presence of groundwater with a total dissolved solid (TDS) concentration of approximately 10,000 mg/L at a depth of 600 m, 16 km west of Biscayne Bay (Williams and Kuniansky 2016). If the TDS concentration is 10,000 mg/L, the chloride concentration must be less than 10,000 mg/L. Thus, the deep water west of Biscayne Bay cannot, on its own, account for chloride concentrations that exceed 18,000 mg/L in the saltwater wedge, and the mechanism by which chloride concentrations might increase along regional-scale flow paths suggested by Weyer (2018) remains unexplained. Furthermore, high-salinity groundwater in areas that might potentially recharge the coastal aquifers of Biscayne Bay are generally associated with mapped low-permeability units. Williams and Kuniansky (2016) find that brackish-to-saline groundwater that lies near the base of the aquifer system in the vicinity of the central part of the Georgia-Florida state line is disconnected and is probably "trapped connate water in fine-grained carbonate rocks near the base of the system isolated from higher permeability rocks above." Such connate water, even if it flowed on a regional scale, would represent a limited, and therefore unlikely, source of sustained saltwater flow to the Biscayne aquifer adjacent to Biscayne Bay. In any case, Weyer's (2018) explanation is, at best, incomplete, and is not the simplest one that fits the available observations. The close similarity between the groundwater chloride concentration and that of seawater, together with the flow pattern inferred by Cooper et al. (1964) from their head and tracer observations, indicates that lateral seawater intrusion from Biscayne Bay is likely the most straightforward explanation. Weyer (2018) suggests groundwater flow can be called "density-driven" only under "hydrostatic" conditions, that is, when a denser or less-dense fluid is submerged in an otherwise hydrostatic host fluid. Furthermore, Weyer (2018) proposes that when forces resulting from the pressure gradient in the host fluid do not exactly counteract gravity, conditions are "hydrodynamic," and the flow must properly be called a "gravitational groundwater flow system." This would seem like a largely semantic distinction but for the additional assertion by Weyer (2018) that in the hydrodynamic case the "buoyancy force" is not directed vertically, as generally understood, but rather is "in the direction of the pressure potential gradient of the host fluid," that is, according to Weyer (2018), buoyancy always acts in the direction of the pressure gradient.

Buoyancy acts along the direction of gravity
Before discussing Weyer's (2018) interpretation of buoyancy, it will be helpful to clarify the relevant terminology. Archimedes' Principle states that, "a body is buoyed up by a force equal to the weight of the displaced fluid" (for example, Whitaker 1968, p. 53), that is, the buoyancy force is equal to the weight of (the force of gravity on) the volume of fluid displaced by the submerged body. For a fully submerged body, the "apparent weight" of the body is the actual weight of the body less the buoyancy force. A body submerged in a fluid of identical density has no apparent weight and remains suspended in the fluid because the buoyancy force on the body equals its weight. A body of greater or lesser density than the surrounding fluid has a positive or negative apparent weight and tends to sink or rise, respectively, relative to the fluid. Thus, the apparent weight is the change in force on the body that arises by virtue of its density contrast with the surrounding fluid. Although Weyer (2018) uses the term "buoyancy force," it is clearly the change in the driving force caused by the density contrast, and not the buoyancy force itself, that he is discussing. We will call the change in the driving force that arises from the density contrast the "buoyant driving force." Weyer's (2018) misinterpretation of buoyancy apparently stems from Fig. 16 of Hubbert (1953) and the accompanying analysis of the driving forces that act on fluids of different density. In Hubbert's (1953) Eq. 39 and a subsequent equation that is not numbered, the driving forces for oil and freshwater are expressed as and where E o and E w are the vectors of force per unit mass on oil and freshwater, respectively; o and w are the densities of oil and freshwater, respectively; g is the gravitational acceleration vector; P is pressure; and ∇P is the pressure gradient. (The driving force per unit mass, i.e., the negative potential gradient, is correctly depicted in Weyer's Fig. 3 as being the vector sum of contributions from gravity and the pressure gradient, as in Eqs. (1) and (2) above.) The difference between the force per unit mass on oil and the force per unit mass on freshwater is then According to (3), the vector E o − E w is in the direction of the pressure gradient. Note that (1)-(3) can be adapted to describe any pair of fluids, for example, freshwater and seawater, by substituting in the corresponding densities. Weyer's (2018) Fig. 4, which is a modification of Hubbert's (1953) Fig. 16, correctly illustrates the relation between forces per unit mass ("resultant forces") on fluids of various densities: freshwater, "water of ocean-type salinity," saturated brine, oil, and gas. All of the resultantforce vectors share a common origin and terminate along a line that is parallel to the pressure gradient. (The resultantforce vector for gas is not depicted as terminating on that line because of the scale of the figure.) Thus, following the rules of vector addition, the difference between any pair of resultant-force vectors is parallel to the pressure gradient, in accordance with (3). It follows that if one were to take the difference in force per unit mass between a given fluid and freshwater as the buoyant driving force, one would conclude, as Weyer (2018) does, that buoyancy acts along the direction of the pressure gradient. That would be incorrect, however, because it would not account properly for the change in potential energy as a fluid of greater or lesser density moves relative to freshwater.
Consider a droplet of oil embedded in a flow of freshwater. The total change in potential energy resulting from an incremental displacement of the oil droplet relative to the surrounding water is the sum of the individual changes in potential energy of the oil and the water. If a volume V o of oil undergoes displacement dx (a vector displacement in any direction), an equal volume V o of water undergoes displacement -dx as it flows around the oil to fill the void left by the oil, and the total change in potential energy is Substitution of (1) and (2) into (4) and recognition that the ∇P terms cancel then gives This total change in potential energy is due to the work done by the buoyant driving force F bd through displacement dx: Comparison of (5) with (6) then implies that the buoyant driving force is If the magnitude of g is denoted by g , the magnitude of the buoyant driving force is F bd = V o o − w g, which is Archimedes' Principle for a fully submerged body: the apparent weight of the body (the buoyant driving force) is the weight of the object ( V o o g ) less the weight of the displaced water (the buoyancy force, V o w g) . Note that our analysis does not presuppose Archimedes' Principle; rather, the principle follows from a correct analysis of forces. Furthermore, (7) correctly states that the buoyant driving force always acts along the direction of gravity, whether the water is in a "hydrostatic" or a "hydrodynamic" state.
If one were to assume, instead, as Weyer (2018) apparently does, that the buoyant driving force per unit mass of oil is given by (3), the buoyant driving force on a mass M o of oil would be given by which incorrectly states that the buoyant driving force always acts along the direction of the pressure gradient. Substitution of (8) into (6) would then give which incorrectly implies that the total change in potential energy is the sum of changes in potential energy of the oil and an equal mass M o (rather than an equal volume V o ) of water. Note that the incorrect expression (8) agrees with the correct expression (7) only if the fluid densities are equal, o = w , or the pressure gradient is hydrostatic with respect to water, ∇P = w g , or both.

Existing density-dependent groundwater simulation codes correctly use Hubbert's force potential
When force is expressed as the gradient of a potential, that potential is called a "force potential," and when velocity is expressed as the gradient of a potential, that potential is called a "velocity potential." Weyer (2018) implies that hydrologists commonly use velocity potentials when they should be using force potentials to model density-dependent groundwater flow. Although Weyer (2018) only briefly mentions the issue and does not tie it into his arguments against the concept of a saltwater wedge, he implies that improper use of velocity potentials is related to assuming an incompressible fluid, and that it is causing widespread misrepresentation of density-dependent flow. Weyer (2018) also calls for replacement of existing groundwater modeling codes with new codes that adhere to his interpretation of buoyancy. SUTRA Provost 2002, version of September 22, 2010) and SEAWAT (Langevin et al. 2008) are two existing codes mentioned in Weyer's (2018) introduction.
In fact, the forms of Darcy's Law used in SUTRA and SEAWAT are consistent with Hubbert's potential (Hubbert 1940), which Weyer (2018) cites as being the appropriate "gravitationally driven force potential." If kinetic energy is negligible, Hubbert's potential, Φ , for a compressible fluid is (Hubbert 1940, p. 802, Eq. 30 where z is elevation, is the pressure-dependent density of the fluid whose potential is being evaluated, and P 0 is the reference pressure at which the potential is zero. (Here, Hubbert's Eq. 30 has been rendered in slightly more rigorous notation by introducing a dummy variable of integration, , which makes the differentiation with respect to P in the next equation easier to follow.) For constant g , the potential gradient is then According to the expression on the far right of (11), the potential gradient is the vector sum of contributions from the pressure gradient and the gravitational acceleration vector, as depicted in Weyer's (2018) Fig. 3. Now, in terms of the potential gradient, the volumetric groundwater flux, q , through a three-dimensional, anisotropic porous medium is a generalization of the one-dimensional flux in Eq. 80 of Hubbert (1940, p. 819): where k is the permeability tensor and is viscosity. Substitution of (11) into (12) then gives or, for velocity, v, where is porosity. Note that although (14) gives an expression for groundwater velocity, it does not involve a "velocity potential." In fact, Eq. (14), which is derived directly from Hubbert's force potential, holds whether the fluid is compressible or incompressible and is precisely the form of Darcy's Law used in SUTRA for fully saturated groundwater flow (see Provost 2002, version of September 22, 2010, p. 19, Eq. 2.19a with relative permeability k r and saturation S w both set to 1). The forms of Darcy's Law used in SEAWAT and another popular code, FEFLOW (Diersch 2014), are ultimately expressed in terms of equivalent head based on a reference density (for example, equivalent (12) q = − k∇Φ, freshwater head), but they also are based on and mathematically equivalent to (13) for fully saturated, single-phase flow (see Guo and Langevin 2002, p. 11, Eqs. 17-19, which are expressed in a form that assumes the principal directions of the permeability tensor are aligned with the orthogonal coordinate system; and Diersch 2014, p. 122, Eq. 3.258 for a single phase with relative permeability k r set to 1) and are, therefore, also consistent with Hubbert's force potential. Sound theoretical treatment of buoyancy in existing codes such as SEAWAT, FEFLOW, SUTRA, and SUTRA-MS (Hughes and Sanford 2004) enables them to simulate successfully a variety of density-dependent groundwaterflow phenomena observed in laboratory experiments and in the field. Benchmark comparisons between experimental and numerical results include systems that exhibit "upconing" of saltwater induced by injection and withdrawal of overlying freshwater (Oswald and Kinzelbach 2004), formation of saltwater wedges Chang and Clement 2012;Badaruddin et al. 2015) and freshwater lenses (Stoeckl et al. 2016), and saltwater fingering (Post and Simmons 2010) and double-diffusive convection (Hughes et al. 2005;Dausman et al. 2010) driven solely by density differences. Such comparisons provide compelling evidence for the suitability of existing codes for modeling densitydependent flow. As far as we know, there is currently no simulation code based on Weyer's (2018) theory of buoyancy, and, therefore, no demonstration that such a code could successfully represent density-dependent groundwater flow. In light of the problem with Weyer's (2018) theory discussed above, we do not expect that such a demonstration is possible, except perhaps for flow systems that are relatively insensitive to buoyancy effects.