Ghislain de Marsily is best known for his contributions to the identification of aquifer heterogeneity [Ackerer et al. 2023], his book [de Marsily 1986], or his recent works on the global relevance of water [e.g., de Marsily 2009], to list a few of his contributions. Less known is his work on seawater intrusion and specifically on the puzzle of coastal springs that discharge water containing a large fraction of seawater above sea level [Arfib and de Marsily 2004; Arfib et al. 2007], which motivates this paper.
Advances in science have overturned the ancient Greeks belief that seawater, consisting of water and earth (two of the five basic elements), flowed inland where its earthen component slowly returned to the ground (its natural state) giving rise to fresh groundwater. It is believed that this concept originated in the Island of Cephallonia (Greece) as a result of the simultaneous observation of a continuous flow of seawater into a sink on the west side of the island, and frequent clouds over the island’s mountains [Stringfield and LeGrand 1969]. The fact that seawater could be seen entering the ground, together with the observation that water got fresher (and lighter) uphill, suggested that by the time it reached the mountain top, it had turned so light as to become air (clouds), “clearly” showing that the hydrological cycle was just the opposite of what we now know. The particular phenomenon of seawater inflow in Cephallonia triggered a scientific controversy that would last for many centuries. A number of imaginative but often physically unacceptable hypotheses had been advanced to account for the phenomenon [Fouqué 1867; Wiebel 1873, in Breznik 1973; Crosby and Crosby 1896 and Fuller 1907] until Maurin and Zoetl  demonstrated with a colouring test that the inflowing seawater returned at brackish karst springs located on the opposite side of the island at a distance of 15 km.
Brackish springs are relatively frequent in coastal carbonate formations and their existence has been extensively reported, especially in the Mediterranean area [Payne et al. 1978, reviews of Breznik 1973; Leontiadis et al. 1988; Fleury et al. 2007; Acikel and Ekmekci 2018]. In fact, more than 300 springs have been identified in the coast of former Yugoslavia alone [Bonacci and Roje-Bonacci 1997]. They essentially consist of inland or submarine karst outlets discharging waters with flow-dependent salinity. They may discharge high flow rates with significant salinities (presumably coming from the sea) at elevations of several meters above sea level, which reveals the especial complexity of their hydraulic and salinization mechanisms (see, e.g., records of rainfall, flow rate, salinity, and temperature at S’Almadrava spring, which discharges at 8 m.a.s.l. in Majorca, Spain, Figure 1a). Although these phenomena have been studied for many years, controversy persists.
The occurrence of a well-developed deep karst system in a seawater intrusion zone appears to be the key factor in the formation of brackish springs. This karst could be formed by the combination of tectonic activity and variations in the sea level over geologic time scales [Fleury et al. 2007]. Moreover, mixing of fresh and seawater may lead to carbonate dissolution [Rezaei et al. 2005], further broadening the conduit section at depth and modifying the seawater intrusion pattern. Under these conditions, a plausible scenario consists of a karst conduit open directly to the sea, allowing mixing at a deep conduit branching [Breznik 1973]. The conduit branching is the junction of three conduits: one from inland with freshwater, another one connected to the sea, and the third one with brackish water leading to the spring mouth (Figure 1b). Such conduit branching has never been directly observed, and the uncertainties associated with its geometry and dimensions have led authors to propose different physical mechanisms to explain salinization at heights above sea level. Proposed salinization mechanisms can be divided into two groups [Stringfield and LeGrand 1969]: salinization due to hydrodynamic effects and salinization due to the greater density of seawater. The first group involves physical conditions of flow in tubular openings or other channels which serve as Venturi tubes. This requires a narrowing of the freshwater conduit at the intersection with the seawater one, allowing the freshwater flowing through to cause a depression, of the order of v2∕2g that sucks in seawater. This mechanism has been proposed for springs showing a rising or steady salinity with increasing discharge flow, such as the Slanac spring (Croatia) [Breznik 1973] and the Makaria spring (Greece) [Maramathas et al. 2006]. However, considering the salinity and elevation of the springs, these situations would require very high groundwater velocities and particular conduit networking morphologies that would be difficult to find in nature or to maintain without erosion, not to mention that mixing fresh and seawater leads to aggressive conditions [i.e., limestone dissolution, Sanford and Konikov 1989; Rezaei et al. 2005; Sanz et al. 2011]. The Venturi effect was also conjectured by McCullough and Land  to explain the inflow of seawater into submarine karst conduits in response to tides and waves. However, the inflow is not in phase with sea level fluctuations, which makes the Venturi effect unlikely. Instead, a hydromechanical effect [Guarracino et al. 2012] explains why the maximum flow rate coincides with the mid-point of the rising tide.
The second group of proposed mechanisms assumes that seawater is sucked at the conduit branching when the hydrostatic pressure in the upwards brackish water conduit is less than that in the seawater tube (that is, when the weight of the brackish mixture up to the spring mouth is less than the weight of seawater up to the sea level). This situation becomes possible if the depth of the branching point is large enough to balance the elevation of the spring mouth with the reduced density of brackish water as compared to seawater. The required depth of the branching point is of several hundred meters, which adds mystery to these peculiar springs. Bakalowicz  attributes deep karst conduits to the Messinian salinity crisis some 5.6 million years ago, when the Mediterranean Sea became isolated from the Atlantic Ocean. This caused the sea-level to drop by at least 1500 m [Ryan 1976], thus lowering the base discharge level and the elevation of karst outlets at large depths below current sea level. In fact, not only the low base level, but also the increased salinity contribute to the development of deep karst conduits. Furthermore, as discussed above, mixing of fresh and salt water leads to calcite subsaturation even when both end members are in equilibrium with calcite and the effect increases with the salinity contrast. Enhanced dissolution coupled to the large drops in sea level elevation during the Messinian salinity crisis may explain why these springs are well documented in Mediterranean shores, but rare elsewhere. Kohout [1965, in Kohout et al. 1977] describe saline springs in Florida. Their observations have been attributed to a thermal convection effect (labelled “Kohout effect”): cold (i.e., dense) seawater enters the aquifer, where it warms up, thus becoming lighter, so that it can rise to a significant height. The relevance of this effect increases with the temperature contrast between the cold and hot waters and with the depth of the cold sea. Both factors are large in Florida [Hughes et al. 2007], but relatively small in the Mediterranean Sea, where we have estimated that the Kohout effect could account for only some 25 cm rise above sea level.
A singular feature of Mediterranean brackish springs is the variation salinity with spring discharge, which has been attributed to variations of freshwater pressure with changes in freshwater flow from the aquifer [Gjurasin 1943, and Kuscer 1950, in Breznik 1973]. This mechanism has been widely accepted for springs that display an inverse relation of flow discharge with respect to salinity, i.e., discharging freshwater during high flows and, below a certain flow-limiting value, increasingly saline waters with decreasing discharge. Examples include Almyros of Heraklion spring (Crete, Greece) [Breznik 1973] Blaz spring (Croatia) [Bonacci and Roje-Bonacci 1997], Pantan [Breznik 1973] and S’Almadrava spring [Sanz et al. 2003].
In some cases, geological constraints suggest that seawater contamination of freshwater flowing through a karst conduit would not be necessarily related to a conduit network connected directly to the sea, but rather to a diffusive seawater intrusion from the porous matrix surrounding the conduit [Arfib and de Marsily 2004]. This situation may occur when a freshwater conduit crosses a saline-intruded fissured matrix zone. Arfib and Ganoulis  performed laboratory experiments demonstrating that a considerable mass exchange may happen under these conditions. Vertical conduits, such as the spring conduit in Figure 1, connected to deep horizontal conduits, such as the fresh and sea water conduits, have been described by Lofi et al. , who link these conduits to karstic activity during the Messinian crisis. While the model of Figure 1 makes geologic sense, it is hard to ascertain whether the sea water conduit is open to the sea or connected through a permeable porous medium. We conjecture that some insight into the most dominant mechanism (open conduit or diffusive connection) in a specific case can be gained from modeling.
There have been attempts to apply numerical models to brackish springs. They have been motivated by the fact that the discharge of those springs represents a precious resource in areas with limited water resources. Modeling proved to be an important tool to test different options of spring development proposed in the last decades [Breznik 1973, 1988; Leontiadis et al. 1988; Bonacci and Roje-Bonacci 1997; Cardoso 1997; Sanz et al. 2003]. Some authors have employed nonlinear analysis or inverse modeling to calculate unit hydrographs and impulse responses from rainfall data in the recharge area [Lambrakis et al. 2000; Pinault et al. 2004]. This type of models, extensively used for freshwater karst springs, use box reservoirs to represent the relationship between an input and an output signal. However, these methods do not consider the physical processes and mechanisms controlling the spring functioning. Maramathas et al.  adopted a different approach based on the mass and energy balance on a hydrodynamic analog, which included three reservoirs flowing from tubes lying adjacent to the spring. Two reservoirs represent the karstic system (two karst subsystems with different depletion period), and the third one emulates the sea. This model assumes the existence of a conduit branching with a conduit open directly to the sea (although this is not simulated), and computes discharge and chloride concentration of the spring using rainfall data as model input. Although the model was successfully applied to the Almyros spring of Heraklion, mixing at the conduit branching and the variable-density turbulent flow in conduits were not considered. In fact, the equations governing variable-density turbulent flow have never been addressed before. We believe that these equations are crucial for a full understanding of the physics of spring functioning. In contrast to this work, Arfib and de Marsily  applied a different conceptual model to the same spring. They assumed the existence of a single freshwater conduit surrounded by a saline porous matrix, where salinization of freshwater in the conduit is a consequence of saline flux from the matrix. This approach considered constant-density turbulent flow in the conduits and mass exchange ratio driven by the head difference between the conduit and the matrix. Surprisingly, these two studies at Almyros spring provided similar results even though they were based on rather different conceptual models and approaches. Therefore, a detailed description of these conceptual models for brackish springs and a direct comparison among them should be addressed.
The objectives of this study are (1) to derive the equations governing density-dependent turbulent flow in open conduits, (2) to describe the salinization mechanisms of inland brackish springs presenting a connection with the sea through an open karst conduit or a diffusive seawater connection, and (3) to compare the spring discharge and concentration response for these two salinization mechanisms.
2.1. Conceptual models of brackish springs
Any conceptual model aimed at explaining brackish springs must include a well developed deep karst system and identify the dominant mechanism of seawater contamination, i.e., the way in which seawater intrudes into the aquifer and mixes with freshwater. Here we analyze two conceptual models that synthesize the features of the most plausible natural systems (Figures 1 and 2). Both conceptual models can explain inverse relations of concentration and discharge of inland springs but apply different salinization mechanisms. The first conceptual model, Turbulent–Turbulent (T–T), assumes that groundwater flows only through a network of conduits and that the seawater contamination occurs at a deep branching point through a conduit of undefined shape connected to the sea (Figure 1). In this environment, ground water follows the hydraulic laws for pipes, resulting on laminar or turbulent flow depending on the velocity and properties of the fluid, and the shape and extent of the conduit section [Chadwick and Morfett 1998]. This problem can be simplified into a mass and energy balance at the conduit branching. Mixing of fresh and saltwater occurs when the interface is placed at the conduit branching (Figure 1a). If the freshwater flow is high enough, freshwater will intrude the seawards conduit and create a submarine spring of, eventually, freshwater (Figure 1b). By contrast, if freshwater flow is very low (or eventually zero), seawater will intrude into the freshwater conduit thus increasing the seawater contamination in the aquifer.
The second conceptual model for spring salinization, which we term Turbulent–Porous (T–P), is a simplification of the model described by Arfib and de Marsily . They assumed that the freshwater conduit is embedded in a fissured matrix intruded by seawater. For the sake of simplicity, we treat this fissured matrix as a conduit where water flows according to Darcy’s Law (Figure 2).
Given the complexity of karst systems, the schemes of Figures 1 and 2 should be viewed as modeling simplifications. Combinations of multiple branching points contaminated at different depths are likely to occur in natural systems. Nevertheless, the simplicity of these two models should facilitate understanding of the physics controlling the phenomenon, which is the aim of our work.
2.2. Governing equations
The conceptual models considered in this study consist both of a mixing place (i.e., conduit branching) in which a conduit connected to the sea and another conduit from the aquifer inland, join a third vertical conduit with mixed water leading to the spring mouth. Depending on the conceptual model, the seawater conduit is assumed to be an open conduit (T–T case, Figure 1) or a porous/fissured medium (T–P case, Figure 2). In both cases flow is governed by mass and either momentum or energy conservation. We have found no equations for variable density pipe flow in the literature. Therefore we derived them below.
We consider a conduit with open area A(l) (A = A′𝜙, where A′ is the total area and 𝜙 porosity, in the case of a porous medium), where l is the length along the conduit axis. Fluid mass conservation is expressed as:
Momentum conservation can be expressed in lagrangian coordinates, in which case it expresses Newton’s second law, or in eulerian coordinates. We adopt the second option, which implies that the variation of momentum (∂𝜌Av∕∂t) is equal to the inflow minus outflow of momentum per unit length of conduit ( − ∂𝜌Av2∕∂l) plus the forces acting on the fluid (pressure, weight, and conduit resistance), which, expressed per unit length equals − ∂PA∕∂l −𝜌 Ag cos𝜃−fp, where cos𝜃 = ∂z∕∂l, z is elevation, and fp represents the component along l of the forces exerted by the conduit walls over the fluid. That is,
This equation can be used, together with Equation (1), for solving the problem. It can also be written in terms of energy conservation, which may be more intuitive. To this end, we expand the time derivative in Equation (2), use Equation (1) to eliminate the resulting ∂𝜌A∕∂t and perform some minor algebraic manipulations to obtain
This equation expresses energy conservation explicitly, but it is still not easy to apply because of the pressure forces exerted by the conduit walls. The normal component of these cancels if A is assumed constant. Moreover, measurements are rarely more frequent than hourly so that pressure waves can be ignored and the fluid assumed to be incompressible. Therefore, the flow rate (Q = Av) is constant along the conduit. With these simplifications, Equation (4) becomes
Equation (5) resembles Bernouilli’s equation except for the fact that 𝜌 varies in space and time in response to variations in the proportion of seawater. As a result, 𝜌 cannot be factored out and energy conservation cannot be written in terms of head (energy per unit weight of fluid), but in terms of energy per unit volume. This can be viewed as e =𝜌 gh. Another difficulty is that e is no longer a potential. That is, flow is no longer given solely by differences in local e, but it depends also on the values of 𝜌 along the path. To illustrate this result, assume steady-state and integrate Equation (5) along the conduit, between points 1 and 2, which yields:
The expression of depends on whether the conduit is open (turbulent flow) or porous (Darcy flow). In the case of the open conduit, we have used Manning’s equation:
The fact that 𝜌 varies in space and time forces us to solve the salt mass conservation equation, which we have written (neglecting dispersion) as:
Solution of these equations requires specifying boundary conditions and continuity at inner nodes (branching point). Pressure is specified as equal to atmospheric pressure (zero relative pressure) at the spring mouth. Energy (per unit volume), , is specified at the seawater entry point, which can be derived from
3.1. TURBOCODE solver
The above equations are solved using an iterative algorithm (called TURBOCODE) programmed in FORTRAN. Basically, input data include the description of the system (elevation, area, etc.) and the freshwater flow rate. Output consists of the seawater flow rate, Qs, and derived variables such as mixed water flow rate and concentration (Qm amd Cm, respectively). The procedure is as follows:
- Initialization. Read Qf(t), zspring, zs, zB, Lm, Ls, ns (or ks), nm (or km), As, Am. Time loop. For each time step perform the following operations.
- Assume initial value for Qs [Qs(t) = Qs(t − 1); Qs(0) = Qf(0)].
- Using Qf and Equations (13), (14) and (11) obtain Qm, cm and 𝜌m, respectively.
- Solve transport using Equation (10) by means of a particle tracking method. Once the spatial distribution of concentrations is known, compute . Repeat for the seawater and mixed water conduits.
- Compute HBm and HBs using Equations (16) and (17), or Equation (18).
- If HBm ≈ HBs, go to next time step. If HBm > HBs, reduce Qs (otherwise, increase Qs), and go to step 3.
3.2. Model setting
Simulations are first performed under steady-flow conditions of constant input freshwater flow to facilitate the understanding of the physics of the problem without memory effects. Transient simulations with a time-dependent freshwater input flow are used later to reproduce the actual functioning of brackish springs. The parameter values used in the simulations are listed in Table 1. They do not respond to any particular brackish spring but are selected partly from S’Almadrava spring (Mallorca, Spain) and partially from literature values of other springs of this type. The elevation of the spring mouth is +8 m.a.s.l. for all the simulations.
Parameter values used in the steady-flow and transient simulations, for Turbulent–Porous (T–P) and Turbulent–Turbulent (T–T) conceptual models
|T–P model||T–T model||T–T model|
|ns (m−1∕3⋅s−1)||-||1.5 × 10−2||1.5 × 10−2|
|nm (m−1∕3⋅s−1)||1.5 × 10−2||1.5 × 10−2||1.5 × 10−2|
See text or Nomenclature list for definition of parameter symbols.
4. Results and discussion
4.1. An energy balance analysis of Freshwater-Seawater mixing at the conduit branching
For a given freshwater flow rate entering the conduit branching, mixing can be viewed as resulting from a mass and energy (or momentum) balance problem. To understand this balance we first calculate the energy at the conduit branching measured from the seawater side (HBs) as a function of seawater flow (Qs), after losses along the seawater conduit (Equations (17) or (18)). This energy is compared to HBm (Equation (16)), energy needed to bring the resulting mixed flow to the spring mouth (that is, the energy at the conduit branching measured from the mixed water side). The problem is solved when both energies become identical (Equation (15)). In fact, this is the way TURBOCODE operates. We consider steady-flow conditions, i.e., a constant density in every conduit, for discussion simplicity.
Results for the T–P conceptual model are shown in Figures 4a and 4b, for two series of small and large Qf values, respectively. HBs decreases linearly with seawater flow rate because energy loss at the seawater conduit obeys Darcy’s law (recall that energy at the sea intake is constant and equal to the hydraulic pressure of seawater). This linear relationship is specific to the seawater conduit and does not depend on how much freshwater is flowing into the system. By contrast, the energy necessary to push a column of mixed water towards the spring mouth, HBm, depends on both freshwater (Qf) and seawater (Qs) flow rates at the conduit branching. For any given Qf, this energy increases with Qs because both friction losses and density of the mixture (and thus the weight of the water column) increase with the mixing ratio. For a low Qf (e.g., 0.25 m3/s, Figure 4a) energy losses in the conduit are small and the dominant term controlling the overall energy is the weight of the water column (Equation (5)), which decreases with decreasing Qs (i.e., with the density of the mixture). For slightly larger Qf (e.g., 1.0 m3/s, Figure 4a) we find that less energy is necessary to push up the same Qs (Figure 4a). This reflects again the fact that energy is most sensitive to the density of the mixed water in the vertical column, which decreases with an increasing Qf/Qs ratio. In fact, HBm is virtually insensitive to Qf for null Qs. This can be attributed to the fact that, for small Qf, friction losses in the vertical conduit are also small and the energy required to push the mixed water up depends only on the weight of the water column, which is insensitive to Qf when Qs is near zero (Figure 4a).
The situation changes for large values of Qf (3.0 to 6.5 m3/s, Figure 4b) as the energy loss associated with flow resistance in the upwards conduit become increasingly important. In fact, above a critical Qf (around 2.5 m3/s for this example) friction losses become more important than density effects. As resistance increases quadratically with the flow rate (Equation (7)), so does the energy required to allow the mixed water to ascend. It should be noted that the increase in energy needed to compensate the resistance of the wall of the conduit is not very sensitive to Qs because for such high values of Qf the proportion of Qs in the total spring discharge is minor.
Every point in which two energy lines (HBm and HBs) cross in Figure 4 represents a solution of the problem for steady flow regimes, i.e., the energy equilibrium point at the conduit branching for a specific freshwater–seawater mixing ratio (Equation (15)). In summary, for increasing Qf values, equilibrium occurs with increasing Qs below the freshwater critical value (Figure 4a), but with decreasing Qs above it (Figure 4b). The overall trend is represented as a freshwater–seawater curve in Figure 5a. The freshwater–seawater curve is particular for every brackish spring and it is representative of the dimensions of the karst system and the physics of seawater contamination at depth. This critical Qf value (around 2.5 m3/s for this example, Figure 5a) separates the flow conditions in which the system is primarily controlled either by the weight of the water column or by friction losses. Note that for very large Qf values (higher than 6.7 m3/s, in this example) the energy equilibrium at the conduit branching is reached for energies higher than the hydraulic pressure of seawater at the conduit outlet to the sea, i.e., for negative Qs. When this happens, freshwater flows through the conduit branching towards the sea and a submarine spring is created. This situation is addressed in detail in Section 4.3.
Figure 5b displays the relationship of salt mass fraction and the spring discharge. The relation reproduces the freshwater–seawater curve discussed above with a decrease on the salinity of the discharge with the flow rate.
For the T–T conceptual model, the behaviour of the mixed water conduit is identical to the T–P case because differences are restricted to the seawater conduit. Still, the T–T case also display two ranges of Qf depending on whether the solution is controlled by the weight of the mixed water column, or by the energy loss in the conduit (Figure 5a). Differences arise in terms of the energy necessary for a particular seawater flow rate to occur (HBs, Equation (17)). The energy loss in the seawater conduit is no longer linear, but quadratic with Qs, and so does the decrease of HBs. As a consequence, the energy equilibrium at the conduit branching is reached for different freshwater–seawater mixing ratios on each conceptual model. Figure 6 illustrates these differences for a low and a high Qf in our examples (0.25 and 4.5 m3/s). Because the examples chosen to illustrate the two conceptual models are designed to hold the same freshwater–seawater ratio at the Qf critical value, any equilibrium point for lower or higher freshwater flow rate yields a lower freshwater–seawater ratio at the conduit branching for the T–P conceptual model than for the T–T case (Figure 5a).
The relation of salt mass fraction and spring discharge in the T–T model differs from that obtained for the T–P case near the extreme values of spring discharge (Figure 5b). Thus, the salinity of the spring for lower and higher flow rates on the T–T model increases slightly with respect to that of the T–P case to finally end at the same point. The relationship of the curve at extreme values of spring discharge may therefore be used as a distinctive feature of T–P or T–T conceptual models. This issue is discussed in detail later.
4.2. Solution for the first extreme point (zero seawater and freshwater flow rates)
We discuss here the two points with zero seawater flow in the freshwater–seawater curve. Figure 5a shows that Qs tends to zero when Qf approaches zero and for a particular high value of Qf. Figure 5a shows that when the freshwater inflow from the aquifer decreases to zero, the seawater intrusion at the conduit branching also decreases but in a lesser degree. As a consequence, the density of the mixed water flowing up towards the spring mouth increases as does the weight of the water column connected to the spring mouth. Eventually, energy loss and kinetic terms become negligible. Ignoring all velocity dependent terms in Equations (16) and (17) or (18) yields
On the contrary, for the T–P case the concentration at the spring mouth decreases linearly with Q for small Q values (Figure 5b). In fact, equating HBm (Equation (16)) and HBs (Equation (18)) and taking derivatives yields
That the behaviour near the origin does not depend much on hydraulic parameters is also shown in the freshwater vs. saltwater discharge curve. Taking derivatives of the mass balance equations (13) and (14) with respect to the Qf yields
In summary, the behaviour of the Qs–Qf and cm–Q curves at the origin is informative of the depth of the branching point, which can be derived either from cm or from the slope of the Qs–Qf curve. It is also informative of the type of seawater conduit: turbulent if the slope of the cm–Q curve is zero, or darcian, otherwise. In the latter case, the slope of the cm–Q curve at the origin allows one to derive the hydraulic resistance of the seawater conduit.
For our example, representing a brackish spring with a spring mouth at +8 m.a.s.l. and a conduit branching at −540 m.a.s.l. (Table 1), we obtain a maximum potential salt mass fraction at the spring mouth of 1.4 × 10−2 (equivalent to 38% of seawater) (Equation (20)). The grey area in Figure 7 illustrates all freshwater–seawater combinations representing a salt mass fraction over the maximum potential calculated, i.e., combinations out of the calculated potential range. In other words, it shows freshwater–seawater mixing ratios that would never be measured at the spring discharge of our modeled brackish spring. As expected, for either conceptual model considered, the solution fits the maximum potential freshwater–seawater ratio for very low Qf values. As the Qf increases the solution separates from that potential ratio given that the energy loss becomes more significant. It should be pointed out that the dependence of Qs on Qf is more linear in the T–T case than in the T–P case. This is attributed to the fact that HBs remains essentially constant for low Qs (Figure 6). This implies that salinity at the spring mouth will tend to be constant value for a range of very low flow rates (say, below 0.2 m3/s in our example) in the T–T case, while it will grow steadily to the same value in the T–P case. It relevant to point out that this quasi linear behaviour does not require a Venturi effect, as suggested by Maramathas et al.  for the Makaria spring.
4.3. Second extreme point (Qs = 0,Qf > 0) and submarine freshwater spring
The second situation in which we have zero Qs occurs for a particular high value of Qf (6.7 m3/s in our example, Figure 5a). As discussed before, the energy loss at the mixed conduit increases quadratically with the mixed water flow, so that the seawater flow reduces to minimize the concentration of the mixed water, and therefore the weight of the mixed water column. This explanation implies that the resulting flow rate is independent of the seawater conduit behaviour. In fact, using again the energy continuity condition (Equation (15) with (16b) for HBm and either (17) or (18) for HBs) together with Qs = 0, yields
Simple inspection of this equation shows that Qf at this point solely depends of the elevation of spring and branching point and turbulent losses in the mixed conduit. Energy continuity, together with fluid and solute mass conservation can be used to compute the derivatives of Cm with respect to Qm and Qs with respect to Qf. The derivation is a bit tedious and the results less informative than at the origin. The only relevant result is that the slope of both curves is far larger in the T–T case, than in the T–P case. This finding is related to the fact that head losses in the seawater conduit are proportional to in the T–T case, so that they are negligible for small Qs. They are proportional to Qs in the T–P case, so that they are sizeable even for small Qs.
For even higher Qf, the functioning of the spring is inverted. The energy necessary at the conduit branching to exceed the energy loss in the mixed water conduit is too high to be maintained. As a consequence, freshwater starts to flow towards the sea reducing the mixed water flow upwards. The freshwater intruding into the conduit connected to the sea promotes a negative Qs (i.e., towards the sea) and a submarine spring is created. Concentration at the spring mouth becomes pure freshwater while the concentration in the submarine spring mouth is initially that of pure seawater. When that happens the functioning of the system changes radically and becomes dependent on the magnitude of the freshwater intrusion into the conduit connected to the sea. In fact, the magnitude of the freshwater intrusion will determine the salinity of the submarine spring, which may even become pure freshwater in an extreme situation. It is important to note that the minimum Qf value for which a negative Qs is produced is essential for a complete characterization of the spring functioning since it defines the flow ranges for which the system behaves as an inland or a submarine spring. It is also important to note that although any freshwater–seawater curve calculated for a particular brackish spring has a theoretical Qf before the Qs is inverted, the actual chances for that high Qf to occur in the aquifer may be very low and that situation could never happen in reality. Measuring pure freshwater at the spring mouth during high flow discharge periods may be an indicator that freshwater is intruding into the seawater conduit thus forming a submarine spring.
To obtain a more realistic solution to the problem involving submarine springs, equations must be solved in transient-flow mode. In order to do that, a time-dependent freshwater inflow function is designed with an increase from 5.0 to 9.0 m3/s followed by a decrease back to 5.0 m3/s, with a variation rate of 0.5 m3/s every 24 h. This freshwater function is introduced in a T–T simulation type with a sea opening at −700 m.a.s.l. (Table 1). Note that for the T–T conceptual model the depth of the sea outlet needs to be defined (Equation (17)), since the average density in the conduit connecting with the sea may no longer be equal to that of the seawater. Results are shown in Figure 8.
Under these flow circumstances, the controlling parameters are the depth of the sea outlet and the average density in the conduit connected to the sea (which depends on the length of this conduit and the magnitude of the freshwater intrusion into this conduit). If the depth of the sea outlet is equal or less than that of the conduit branching, the pressure at the sea outlet is therefore less than that at the conduit branching, and freshwater is able to easily intrude into the conduit. The submarine discharge will become fresh very quickly, only depending on the length of the conduit.
By contrast, if the sea outlet is deeper than the conduit branching (as it happens in our example, Figure 8), the velocity of the freshwater intruding into the seawater conduit will be very small because pressure at the branching point needs to increase to overcome the double effect of having a higher pressure at the sea opening, and the reducing density of the water filling the conduit, initially assumed to be pure seawater. In Figure 8, freshwater flow, Qf, increases from 6.7 to 8.75 m3/s, with a negligible flow in the seawater conduit. As Qf increases further, pressure at the branching becomes sufficient to overcome the reduced density of freshwater filling the sea conduit and the submarine spring discharges pure freshwater (freshwater flow above 8.75 m3/s, Figure 8). When Qf decreases again (e.g., after an important rainfall event) the behavior of the system inverts and the fluid in the conduit becomes positive again (freshwater flow reduction from 8.75 to 6.7 m3/s, Figure 8). Note that for this range of flow rates, although groundwater flows back towards the conduit branching, the mixed water moving upwards to the spring is still pure freshwater. Only when the seawater front reaches the conduit branching (at 6.7 m3/s approx., in our example) will the spring discharge become salinized again.
4.4. Sensitivity analysis
A set of simulations are conducted to illustrate the sensitivity of the spring response with respect to the parameters used in the simulation. First of all, we identify the controlling parameters for each conceptual model on steady-flow conditions from Equations (16) and (17) or (18). Then a perturbation of ±20% is applied separately over the base value of every parameter in order to quantify the relative sensitivity of the solution to that parameter (see parameter unperturbed base values in Table 2). The analysis is performed over steady-flow simulations since it allows a better isolation of the effect that each parameter has on the solution. The sensitivity analysis is discussed separately for the two conceptual models. Results are presented in terms of freshwater–seawater curves and relation of salt mass fraction with spring discharge. These representations are especially useful because they can be compared with the most commonly available field measurements in brackish springs, which are spring flow rates and water concentration with time.
Values of the controlling parameters used in the base simulations for the sensitivity analysis, for Turbulent–Porous (T–P) and Turbulent–Turbulent (T–T) conceptual models. These values are perturbed by ±20% to complete the analysis
|T–P model||T–T model|
|Am (m2)||2.0||Am (m2)||2.0|
|nm||1.5 × 10−2||nm||1.5 × 10−2|
|zB (m)||540.0||zB (m)||540.0|
|Ls∕ks As (m−3)||6.3 × 10+7||(m−2.5)||1.41|
See text or Nomenclature list for definition of parameter symbols.
4.4.1. Turbulent–Porous conceptual model
The solution for a T–P problem depends only on four parameters: Am, nm, Ls∕ks As and ZB. Although this analysis is performed for steady-flow conditions, these four parameters (plus the history of Qf(t)) would control a transient-flow simulation as well. Figure 9 displays the results (Qs vs. Qf, and cm vs. Qm) for the base case and perturbated parameters.
The parameters Am and nm control the resistance of the vertical conduit to the water flow (Equation (16)) and therefore show a similar—but opposite-influence on the spring response (Figure 9a–d). As the energy loss in the vertical conduit is a quadratic function of the water flow, the effect of these parameters is negligible for low Qf or spring discharges (recall that for low Qf the solution was primarily dominated by the weight of the water column, Figure 5a). However, their influence rapidly increases with Qf (Figure 9a,c). As a consequence, the expected spring concentration will increase with decreasing flow resistance, although this effect can be observed only for high spring discharges (Figure 9b,d). Am not only controls the energy loss but also the velocity of the groundwater in the conduit, and therefore for transient-flow situations presents a relatively stronger influence on the spring response than nm.
The solutions display a negative sensitivity with respect to the parameter Ls∕ks As (Figure 9e,f), i.e., an increase in the resistance to flow in the seawater conduit reduces the proportion of seawater entering the conduit branching. Obviously, the solution does not change near the two zero-Qs situations, which is also consistent with the equations governing the functioning mechanism of the spring. As discussed above, when Qf approaches zero, the solution only depends on the depth of the conduit branching (Equation (18)) and, accordingly, the solution shows no sensitivity to any other parameter than zB (Figure 9g,h). Further, for Ls∕ks As the solution shows no sensitivity at the other zero-Qs point either because the energy equilibrium for this situation to occur depends on the depth of the conduit branching and on the resistance in the vertical conduit (Figure 9, Equations (16) and (18)).
Finally, the depth of the conduit branching, zB, exerts an influence on the spring concentration for any discharge (Figure 9h). As discussed above, the proportion of seawater entering the conduit branching increases with zB because the pressure resulting from the weight of the seawater column increases accordingly (Figure 9g). Consequently, the spring salt mass fraction for any spring discharge is higher for deeper conduit branchings (Figure 9h). This parameter influences the energy curves of both the seawater and the mixed water conduit (Equations (16) and (18)). The slope of the freshwater–seawater curve when the freshwater flow approaches zero, is also dependent on the parameter value (recall Equation (19), Figure 9g). Results also show that shallow conduit branchings will facilicitate the formation of submarine springs because the minimum Qf necessary to inverse Qs is smaller.
An overall analysis of Figure 9 suggests that a spring with measures in a broad range of discharges may allow us to characterize: (1) the resistance of the upwards conduit (but not nm and Am Lm separately) from the response at high flow rate; (2) the depth of the branching point (zB) from the concentration at low discharges; and (3) the hydraulic resistance of the seawater conduit from the slope of cm vs. Qm at low spring discharges.
4.4.2. Turbulent–Turbulent conceptual model
The solution for a T–T problem depends on four parameters: Am, nm, and ZB. Figure 10 shows the results of the sensitivity of the solution to these parameters. Results for this case are generally analogous to those of the T–P case and similar discussion could be applied. The most striking feature is that the behaviour at the two zero Qs points is virtually insensitive to flow resistance parameters at the seawater conduit. This is consistent with Equations (21) and (24). The main difference with the T–P case is that results are less sensitive to zB and that the abrupt change in Qs near the critical Qf remains nearly vertical for all parameter combinations (reflecting that energy losses in the seawater conduit are proportional to , so that they drop more significantly to zero when Qs tend to zero than in the T–P case). The other relevant feature is that, for all parameter sets the maximum Qs is flatter in the T–T case than in the T–P case.
In summary, the sensitivity analysis confirms that the main features to distinguish the T–T and T–P models are: (1) near the origin: the slope of the cm vs. Qm curve (near zero in the T–T case), (2) near the critical Qf point, the slope of the Qs vs. Qf curve (nearly vertical in the T–T case), and (3) near the maximum Qs point, the flatness of the Qs vs. Qf curve (inverted parabola and asymmetric in the T–P case and flatter, nearly circumferential in the T–T case).
4.5. Comparison with field data
The results obtained in TURBOCODE simulations are compared with field data observations from three brackish springs: S’Almadrava spring (Mallorca, Spain, Figure 1), Pantan spring (Croatia) and Almyros of Heraklion spring (Crete, Greece). Figure 11 shows the relationship of the spring concentration vs. discharge for S’Almadrava [Sanz et al. 2002], Pantan [Breznik 1973] and Almyros [Monopolis et al. 1995]. Beyond differences in number of observations and in salinity or discharge ranges, all springs display similar patterns with a decreasing salinity for increasing spring discharges (Figure 11). Based on the data available, the authors proposed that the salinization of these springs is due to the greater density of seawater in the T–T scheme, as discussed below. Field data for S’Almadrava covers a wide range on spring discharges (Figure 11a). From these observations, the conceptual model controling the functioning of S’Almadrava spring would be a T–T case because: (1) the relation of salinity and flow rate is far from linearity for medium spring discharges, (2) for very low discharges, the spring concentration seems to reach a constant value (Figures 5b and 10, right hand side plots), and (3) Qs remains quite constant near its maximum.
The available observations for Pantan spring also cover a wide range of spring discharges and the connection to the sea seems to be controlled by a karst conduit (T–T case) (Figure 11b). However, the limited number of data available prevents further discussion. In contrast, many observations are available for Almyros spring (Figure 11c). Data concentrates for medium to high spring discharges, where the spring clearly discharges freshwater. Field measurements seems to form a curve of high linearity. This observation seems more consistent with a T–P case (Figures 5b and 9, right hand side figures). However, the distinction of the conceptual model of Almyros turns out to be difficult due to the lack of data for very low flow discharges.
The fact that field measurements display some spreading with respect to the simulated curves for steady-flows presented in this study could be attributed to the memory effects resulting from the dynamic nature of the systems or be a consequence of the higher complexity of the karst network in the aquifers (e.g., multiple conduit branchings), compared with the simplifications considered in this study.
The freshwater seawater flow rates plots from field measurements in S’Almadrava spring and Almyros spring [Maramathas et al. 2006] are shown in Figure 12. The data from S’Almadrava spring shows a zero-Qs point for Qf approaching zero. The freshwater–seawater ratio then increases linearly with Qf before gradually stabilizing in about 0.15 m3/s of Qs for a wide range of Qf. This pattern agrees with simulations of this study for T–T conceptual model (Figures 5a and 10, left hand side figures). The field data for this spring does not show any significant decrease in Qs for higher Qf, and therefore Qf of at least 4 m3/s are expected to be necessary for the formation of a submarine spring.
Field measurements for Almyros spring form a curve with an abrupt decrease in Qs with Qf until a zero-Qs point is reached for about 5.6 m3/s Qf. Null seawater contribution (that is, pure freshwater discharging at the spring) is maintained for higher Qf values, in line with the formation of a submarine spring at the sea bottom predicted in our simulations (Figure 8). The formation of submarine springs may be reconcilable with either T–T or T–P conceptual models and we would need more geology information to confirm any conceptual model. However, the abruptness of the freshwater–seawater curve when approaching the zero-Qs point indicates that a T–T scheme is more likely at Almyros spring. The different interpretations that may arise from the available observations for Almyros spring point out the difficulty of selecting the conceptual model governing brackish springs. This is especially so when the available observations do not cover a wide range of spring discharges including the two zero-Qs points. In fact, Arfib and de Marsily  and Maramathas et al. , built numerical simulations describing the T–P and T–T conceptual models, respectively, and obtained reasonable results in both cases.
In this study we use hydraulic models to reproduce the physics controlling the salinization of brackish springs. We derived the equations governing variable-density turbulent flow in conduits and solved them numerically. The solution has been found for systems consisting of a deep mixing place in which a freshwater conduit and another conduit connected to the sea join a third conduit with mixed water leading to the spring mouth. Two conceptual models have been tested: the Turbulent–Turbulent case in which the salinization of the conduit branching comes from a karst conduit connected to the sea, and the Turbulent–Porous case in which the seawater intrusion occurs through a porous/fissured matrix (Figures 1 and 2, respectively). In both cases flow is governed by mass and either momentum or energy conservation. The response of the spring concentration to the variation of freshwater flow rate from the aquifer is evaluated in steady-flow conditions in terms of energy balance at the mixing point and freshwater–seawater ratios. For low freshwater flow rates, the spring response is controlled by the weight of the water column flowing from the conduit branching to the spring mouth. As the freshwater flow rate increases, the energy loss increases accordingly and the resistance of the conduit walls to flow becomes the controlling factor. Results are similar for both conceptual models for medium spring discharges, although the relation of spring salinity and discharge is more linear for the Turbulent–Porous case. Additionally, for the Turbulent–Turbulent case, the spring concentration becomes constant when the spring discharge tends to zero. This points out that in this situation the solution depends only on the water column weight in the vertical conduit. In fact, the spring dries up when the excess in elevation of the spring mouth is balanced by the excess of density of seawater with respect to that of the mixed water in the conduit connected to the spring mouth. Further, the depth of the conduit branching can be then approximated from the elevation of the spring mouth and the spring salinity at very low discharges.
The sensitivity of the spring response with respect to the parameters controlling the system is also evaluated in terms of freshwater–seawater ratios and the relation of spring discharge and salt mass fraction. For the Turbulent–Porous case, the simulation is sensitive to Am and nm for high freshwater flow rates, but their effect is opposite: salinity increases for larger Am and smaller nm. This is attributed to the quadratic increase of energy loss in the vertical conduit with the spring discharge. The solution for medium freshwater flow rates depends on Ls∕ks As. A higher parameter value increases the energy loss of the seawater intruding into the aquifer and thus reduces the salinitzation at the conduit branching. Finally, for low freshwater flow rates the solution only depends on the depth of the conduit branching since it controls the limiting water column weight that the system can support before the spring dries up. The sensitivity of the solution in the Turbulent–Turbulent case is very similar to that for the Turbulent–Porous one for the parameters Am, nm and ZB. The resistance to flow in the conduit connected to the sea can be expressed now as ., where an increase of the resistance not only reduces the proportion of the seawater mixing at the branching point but also the salinity discharging from the spring. This effect is only important for medium freshwater flow rates since the solution for high flow rates is dominated by the energy loss in the vertical conduit rather than in the conduit connected to the sea. If the freshwater flow rate overpasses a limiting value, the high energy loss in the vertical conduit may promote a freshwater intrusion into the conduit connected to the sea. This will generate the formation of a submarine spring in the sea floor of variable density depending on the extent of the intrusion and the depth of the conduit outlet in the sea. This situation may be identified by measuring freshwater concentration during stages of high spring discharges in the spring.
The prediction curves of freshwater–seawater and the relationship of spring discharge and salt mass fraction are compared with field data from S’Almadrava spring, Pantan spring and Almyros of Heraklion spring. The simulation results show a good agreement with the field data available, and provide insights to identify the conceptual model governing every particular spring. The analysis highlights the importance of acquiring high frequency field data encompassing the whole range of spring discharges for a sound understanding of the spring functioning, specially at the extreme points of low seawater flow.
The freshwater–seawater curve is specific to every brackish spring and it is representative of the dimensions of the karst system and the salinization mechanism. In the light of our findings, an analysis based on the freshwater–seawater ratios rather than on the relationship of spring concentration and discharge proves to be more suitable for identifying the salinization mechanisms of brackish springs.
|Al||Cross sectional area of a conduit|
|A′||Total area of a conduit|
|l||Length along the conduit axis|
|𝜌||Water density, defined in Equation (11)|
|𝜌0||Density of pure water|
|fP||Shear component of the forces exerted by the conduit walls over the fluid|
|Q||Water flow rate at conduit|
|zspring||Depth at spring mouth|
|zs||Depth at the seawater entry point|
|RH||Hydraulic radius of the conduit (ratio of A to wet perimeter)|
|c||Salt mass fraction|
|H||Energy per unit volume|
|HBf||Energy at the conduit branching measured from the freshwater side|
|HBs||Energy at the conduit branching measured from the seawater side|
|HBm||Energy at the conduit branching measured from the mixed water side|
|f||Subscript for freshwater and freshwater conduit|
|s||Subscript for seawater and seawater conduit|
|m||Subscript for mixed water and mixed water conduit|
|B||Subscript for conduit branching|
|Overbar stands for spatial average|
Conflicts of interest
Authors have no conflict of interest to declare.
Comments by two anonymous reviewers and especially John M. (Jack) Sharp are gratefully acknowledged. Most of this work was funded by the European Union (SALTRANS project, contract no. EVK1-CT-2000-00062) and the Spanish Government (project Medistraes3. grant no. PID2019-110311RB-C21), and IDAEA-CSIC Center of Excellence Severo Ochoa (Grant CEX2018-000794-S). ES is also grateful to a FI doctoral scholarship from the Catalan Government. The authors wish to acknowledge Dr. Guangjing Chen for his help and useful comments on the formulation and code design. Data and documentation related to S’Almadrava spring was gently provided by the Junta d’Aigües de Balears.