## Abstract

Karstic spring discharge is related to the hydraulic head recession through a power function with an exponent <1. In the literature, analytical solutions are available for exponential and non-exponential models based on a set of restrictive physical and mathematical assumptions. The models search for a holistic and deductive solution without basic physical and practical interpretations, simple logical inferences leading to mathematical analytical or empirical formulations. In this paper, an inductive, logical, practical, and instead of holistic modeling, physically plausible piecewise solutions are proposed with detailed inferences and interpretations. In the proposed methodology, the discharge and hydraulic head records are decomposed first into a set of verbal classes and, subsequently, physical meaning for each class is explained leading to simple general but empirical models. For this purpose, Wakula and St. Marks River (Florida) hydrograph records are used for the general solution sinkhole discharge and hydraulic head variations. The solution methodology presented in this paper does not make any distinction between relatively small or large sinkhole heads. The calibration and verification of the methodology is shown with a comparison of the available record values to partial power models. Finally, it is concluded that the proposed methodology is reliable and can be applied to hydraulic head availability with discharge records in any part of the world for karstic aquifer domains.

## INTRODUCTION

Karstic terranes are geological formations where carbon dioxide (CO_{2})-rich waters cause and create cavities, porous parts, fissures, and fractures in the water-bearing layer. Rainfall waters take CO_{2} from the atmosphere, and generate carbonic acid (H_{2}CO_{3}) in the groundwater, which circulates through the fissures, fractures, and cavities and enlarges the karstic voids. Compared to porous and fracture rock aquifer hydrogeological studies, karstic terrane poses difficulties in their flow modeling. They can be visualized as double-porous media for analytical treatments under a set of restrictive assumptions. Geologically, karstic aquifers take place in limestone and dolomite layers, which are soluble media. Ford & Williams (2007) stated that sinkholes and shallow trenches on the karstic domain surface are the main conduits for groundwater recharge. Karstic aquifers, which constitute the only source of water supply in some semi-arid regions, are essential (Pla *et al.* 2016). In many parts of the world (Albania, Austria, India, Bosnia & Herzegovina, China, Slovenia, Turkey, USA, etc.) karstic aquifers provide domestic water supply. Endemic species and rare elements enrich their quality. In the USA almost 25% and in some European countries up to 50% of domestic water is suppled from the karstic aquifers (Elliot 2000; Hartmann *et al.* 2014).

The abstraction of groundwater from karstic aquifers is comparatively easier than from the porous or fracture types because groundwater flows freely within the interconnected network of conduits. However, such free flow invites easy groundwater pollution contamination and movement along long distances in short time durations, which is an unwanted property of karstic aquifers due to their vulnerability to contamination (Field & Nash 1997; Goldscheider *et al.* 2010; Lauber *et al.* 2014). Li *et al.* (2016) discussed the municipal wastes entrance into Wakulla Springs in northwest Florida through Ames Sink, where the leaching of wastes into karstic waters with algae populations, which consume oxygen in the water, caused reduced levels of oxygen in the water, resulting in the deterioration of groundwater quality.

In conduits the groundwater movement takes place rapidly in a turbulent manner, but in the surrounding porous, fissured, and fractured media water flow is laminar. These flow properties give karstic aquifers double-porosity characteristic. Due to the piezometric head difference, an instantaneous response takes place as groundwater recharges through karstic conduits, but the response is rather slow through the fissures, fractures and porous media. Still there is room for karstic aquifer hydrogeological parameter estimations, including the storage and transmissivity values. The complex structure of the subsurface voids (in the form of cavities, fissures, fractures, and porous parts) and the non-linearity of flow laws makes analytical analysis almost unmanageable. For instance, the well-known Darcy Law is not valid in the karstic aquifer parameter quantifications, but non-linear flow laws are valid (Şen 2014, 2019). As mentioned earlier nearly half of domestic water supply is supported from karstic domains in most countries, and, their physical quantifications are indispensable. Any new development concerning these aquifers is therefore welcome in the literature. Since neither porous nor fracture aquifer physical quantification methodologies are valid in the karstic aquifers, it is recommended to monitor their behaviors through a set of main pumping and observation wells or in sinkholes and swallow holes, which are the most efficient way of characterizing their hydraulics for long term planning. Based on the records from these structures, it is possible to explore the holistic behaviors of karstic aquifers as a basis for improved and sustainable management.

At some places it is possible to measure discharges (e.g. from one or more springs) and the hydraulic heads close to the discharge can provide basic information about the physical behavior of the karstic aquifer areal features. The simplest model is the relationship between discharges and water-level variations, such as in the linear reservoir model (Forchheimer 1901; Maillet 1906; Nash 1960). In this manner, it is possible to estimate karstic aquifer water volume availability as proportional to the spring fluxes (Li *et al.* 2016). Of course, spring discharges are not because of point loading, but regional infiltration after each rainfall event. In any mathematical modeling, the infiltration may be thought as the areal average values, and hence, these can be incorporated in a simple analytical model with the assumption of uniform infiltration consideration. In any karstic aquifer there is a set of sinkholes through which groundwater recharge take place. The sinkholes or swallow holes are recharge points and the springs are discharge points in karstic terrane. In general, recharge flow takes place vertically and is converted to horizontal type by the conduits leading to point spring discharges. The holistic mathematical treatment of karstic aquifers is presented through mathematical formulations under a restrictive set of assumptions but one type of model could not be obtained (Li *et al.* 2016). Şen (2018) provided additional information based on their study and showed that rather than a square root exponent a smaller exponent value is more suitable for the karstic aquifers.

There are tremendous research possibilities in the karstic domain and sinkholes as suggested by Li *et al.* (2016), which ignited the research potential of complex karstic aquifers and sinkhole conduits. Their reply (Field *et al.* 2018) to Şen (2018) also clarified any possible misunderstandings between their paper and Şen's (2018) discussion. It is also stated that study of karstic conduits is a challenging task for scientists. Although fluid dynamics should be considered and utilized whenever possible as a useful tool for such studies, basic physical and practical aspects should also provide additional support and useful information for such tools.

It is the main purpose of this paper to provide logical inferences leading to empirical formulations for karstic aquifer spring discharge responses to rainfall occurrences. This methodology simplifies complicated mathematical treatments and provides solution by partial power law relationships about the sinkhole discharge and its hydraulic head variations. The conclusions are compared with the holistic solution results. It is seen that partial verbal empirical approaches provide more scientific information than the existing holistic methodologies. For the application two spring hydrographs are used from the St. Marks Karst Watershed in northwest Florida, which have been already used in several papers in the literature.

## METHODOLOGY

In this paper the sinkhole head, as already used by Li *et al.* (2016), refers to the hydraulic head at the sinkhole and is shown notationally by h_{S}. In their work, the spring location elevation is taken at zero hydraulic head, h = 0, and the sinkhole base piezometric level (hydraulic head) as h_{SB} prior to any rainfall. Spring discharge hydrographs are measured for classical hydrodynamic analysis of the karstic aquifers (Kovács *et al.* 2005; Geyer *et al.* 2008). The rising limb of the hydrograph starts with the precipitation almost immediately. Once the rainfall stops, the peak is reached and then, according to the karstic aquifer physical structure, the recession limb starts to take place. It is possible to make interpretations from the shape of the discharge hydrograph and hydraulic head temporal variations measurements.

Although the groundwater flow continuity equation coupled with a convenient flow law (Forcheimer, Escande, cubic law, etc.) leads to the groundwater balance equation in the form of partial differential equations (which can be solved at a regional scale only by numerical calculations), they need aquifer parameters (storativity and transmissivity), which are difficult to obtain in the karstic aquifer cases. Additional difficulties lie in the spatial configuration of the aquifer, stability analysis satisfaction, tedious calculations, and cost of computation.

_{SP}, is in square root proportionality with the hydraulic head, h

_{S}. A digression is assumed here for more general solutions where instead of square root power (0.5), as was stated by Şen (2018), the exponent,

*n*, in such a relationship should have values 0< n < 1. Hence the basic expression can be written as:where is the hydraulic head and represents the spring flow. It is possible to rewrite this expression by definitions of and and hence:

This last expression cannot be solved by analytical methodology except when *n* = 0.5 for which the solution is presented by Li *et al.* (2016), which will not be repeated here. The right-hand side can be integrated in temporal domain easily, but the left-hand side cannot be integrated by means of the classical methodologies. Although the finite difference and element numerical methods are useful for the spatial integration calculations, they require extensive computation time apart from the stability condition satisfaction.

It is important at this stage that x = 1 leads to an indeterminate state without solution. In Figure 1 various versions of for *n* = 0.1, 0.3, 0.5, 0.7, and 0.9 are presented graphically.

In this paper, instead of holistic mathematical solutions with restrictive assumptions, the available discharge and hydraulic head records are inspected first visually and then, depending on their patterns, the whole record is broken down into a set of few classes, each with distinctive characteristics.

## APPLICATION

The view taken in this paper is that any holistic analytical solution cannot be valid universally for all hydrogeological (karstic aquifers and sinkholes) formations under the light of restrictive assumptions. It may be more illuminating to visualize first some possible partial features in each record and then try to find physically plausible practical solutions for applications. Of course, analytical solutions theoretically based on the mass balance (continuity) and groundwater movement principles with a set of restrictive assumptions are valid approaches leading to general holistic representations. With a set of restrictive assumptions for the whole system, such as homogeneity, isotropy, and uniformity, the partial characteristic details may be important from a physical and practical viewpoint. For instance, Figure 2 (from Li *et al.* (2016)) presents holistic modeling results in the form of exponential and non-exponential (power function) solution, and although the exponential model seems suitable for ‘low’ discharge values, the non-exponential model is suitable for ‘very high’ discharges with no model match in between. The solutions are based on *n* = 0.5 exponent in Equation (1), which may not be valid holistically (Şen 2018). For a better holistic solution, the view taken in this paper is to take either arithmetic or weighted averages of the two holistic models. For this purpose, in the same figure such a model mixture solution is presented with broken red lines. It is obvious that the average model is more representation of the recession discharge parts of records holistically. However, some details are overlooked by the holistic approaches, which consider the aquifer as homogeneous and isotropic media, which is not possible in reality, especially in karstic aquifers.

It is obvious from this figure that none of the models confirm with the general trend of recession. This is because the models try to represent holistically all the recession part in one expression, which cannot be achieved easily. One can visually suggest that either the arithmetic or weighted average mean of the two models provide better results. It is shown with the dotted broken line in Figure 2 and one can appreciate the importance of such an average model, if the problem is to represent the recession limb by means of a single and holistic mathematical approach.

On the other hand, a close look at the discharge record in Figure 2 shows that there are two parts, namely non-linearity at high and medium discharge ranges, but the tail of the record appears in the form of a linear line on the average. This is the main reason why the holistic mathematical approaches cannot be achieved by means of a single mathematical function.

However, one can also appreciate that the sinkhole response to hydraulic head difference is not uniform everywhere. Physically the upper parts of any karstic media have more fractures with porous cavities than deeper positions, where the fine fractures and fissures play a more effective role around sinkholes. On the other hand, the higher the hydraulic head, the faster the water level drop, which is expected to weaken with time and pressure (hydraulic head) drop.

Based on the last two paragraphs, it is easier to appreciate that the recession limb may be thought of as having different parts that are not like each other. For instance, Figure 3 presents some sub-sections of the recession limb with different features. On the left-hand side their descriptions are given as ‘Very high’, ‘High’, ‘Medium’ and ‘Low’ discharge classes.

This figure shows four non-overlapping parts. Each section is shown by straight lines with very different characteristics, which can be summarized as follows:

- 1.
At the ‘V. high’ section the discharge starts to decrease by a certain slope.

- 2.
At the ‘High’ section the slope of the straight line increases, which implies that the discharge change by time is higher than the previous peak adjacent ‘V. high’ class.

- 3.
In the ‘Medium’ section the slope becomes even higher giving the impression that there is a sudden discharge decrease, which may be due to a decrease of karstic aquifer transmissivity reduction, suggesting that the hydraulic head reduction is also sudden, as is obvious in Figure 4.

- 4.
The last ‘Low’ discharge recession limb segment has the smallest slope and the rhythms around the average straight line is meaningful in that the periodic fluctuations become smaller in scale at very low discharge values.

This figure also reflects four distinctive hydraulic head recession classes with different slopes. Each one has a negative slope, but the maximum amount appears in the ‘Medium’ sector. Almost along each sector, there are fluctuations around the straight lines, which indicate the heterogeneity of the Crawfordville (Wakula River) karstic aquifer. At very small scales each fluctuation implies rising part discharges from the nearby fissures, fractures, and/or solution cavities. During the falling sections, the groundwater level in the sinkhole has passed these features and entered comparatively less karstifized locations (Şen 1995; 2018).

Class . | Crawfordville (Wakula River) discharge . | |||||
---|---|---|---|---|---|---|

Start . | End . | Slope, S . | b . | a . | Verification . | |

‘V. High’ | (357, 101.7) | (618, 75.61) | −0.100 | 0.54 | 2,433.80 | 101.70 |

‘High’ | (618, 75.619) | (725, 54.93) | −0.193 | 2.00 | 28,872,581.25 | 74.71 |

‘Medium’ | (725, 54.93) | (784, 34.83) | −0.341 | 5.82 | 2.4371E+18 | 53.82 |

‘Low’ | (784, 34.83) | (1,884, 16.17) | −0.017 | 0.88 | 11,886.23 | 34.83 |

Class . | Crawfordville (Wakula River) hydraulic head . | |||||

Start . | End . | Slope, S . | d . | c . | Verification . | |

‘V. High’ | (208, 2.16) | (549, 2.04) | 0.000 | 0.06 | 2.96 | 2.16 |

‘High’ | (549, 2.04) | (730, 1.83) | −0.001 | 0.38 | 22.60 | 2.04 |

‘Medium’ | (730,1.83) | (806, 1.60) | −0.003 | 1.36 | 13,981.56 | 1.83 |

‘Low’ | (806, 1.60) | (1,555, 1.15) | −0.001 | 0.50 | 46.20 | 1.60 |

Class . | Newport (St. Marks River) discharge . | |||||

Start . | End . | Slope, S . | b . | a . | Verification . | |

‘High’ | (15, 20.84) | (16, 17.1) | −3.740 | 3.06 | 83,819.96 | 20.84 |

‘Medium’ | (16, 17.1) | (20, 14.3) | −0.700 | 0.80 | 157.74 | 17.10 |

‘Low’ | (20, 14.3) | (27, 13.03) | −0.181 | 0.31 | 36.19 | 14.30 |

Class . | Newport (St. Marks River) hydraulic head . | |||||

Start . | End . | Slope, S . | d . | c . | Verification . | |

‘High’ | (15, 2.08) | (16, 2.02) | −0.060 | 0.45 | 7.10 | 2.08 |

‘Medium’ | (16, 2.02) | (20, 1.98) | −0.010 | 0.09 | 2.59 | 2.02 |

‘Low’ | (20, 1.98) | (26, 1.96) | −0.003 | 0.04 | 2.22 | 1.98 |

Class . | Crawfordville (Wakula River) discharge . | |||||
---|---|---|---|---|---|---|

Start . | End . | Slope, S . | b . | a . | Verification . | |

‘V. High’ | (357, 101.7) | (618, 75.61) | −0.100 | 0.54 | 2,433.80 | 101.70 |

‘High’ | (618, 75.619) | (725, 54.93) | −0.193 | 2.00 | 28,872,581.25 | 74.71 |

‘Medium’ | (725, 54.93) | (784, 34.83) | −0.341 | 5.82 | 2.4371E+18 | 53.82 |

‘Low’ | (784, 34.83) | (1,884, 16.17) | −0.017 | 0.88 | 11,886.23 | 34.83 |

Class . | Crawfordville (Wakula River) hydraulic head . | |||||

Start . | End . | Slope, S . | d . | c . | Verification . | |

‘V. High’ | (208, 2.16) | (549, 2.04) | 0.000 | 0.06 | 2.96 | 2.16 |

‘High’ | (549, 2.04) | (730, 1.83) | −0.001 | 0.38 | 22.60 | 2.04 |

‘Medium’ | (730,1.83) | (806, 1.60) | −0.003 | 1.36 | 13,981.56 | 1.83 |

‘Low’ | (806, 1.60) | (1,555, 1.15) | −0.001 | 0.50 | 46.20 | 1.60 |

Class . | Newport (St. Marks River) discharge . | |||||

Start . | End . | Slope, S . | b . | a . | Verification . | |

‘High’ | (15, 20.84) | (16, 17.1) | −3.740 | 3.06 | 83,819.96 | 20.84 |

‘Medium’ | (16, 17.1) | (20, 14.3) | −0.700 | 0.80 | 157.74 | 17.10 |

‘Low’ | (20, 14.3) | (27, 13.03) | −0.181 | 0.31 | 36.19 | 14.30 |

Class . | Newport (St. Marks River) hydraulic head . | |||||

Start . | End . | Slope, S . | d . | c . | Verification . | |

‘High’ | (15, 2.08) | (16, 2.02) | −0.060 | 0.45 | 7.10 | 2.08 |

‘Medium’ | (16, 2.02) | (20, 1.98) | −0.010 | 0.09 | 2.59 | 2.02 |

‘Low’ | (20, 1.98) | (26, 1.96) | −0.003 | 0.04 | 2.22 | 1.98 |

*n*= 0.5 as used in many analytical studies.where f = b/d and e = a/c

^{f}.

The start and end coordinates are given in Table 1 for each straight-line segment. The verification is obtained by calculating the start point values from the empirical expressions. It is obvious that there are relative errors less than ±1%. The application of the empirical model with parameters to the Crawfordville (Wakula River) discharge case yields the model pattern in Figure 3. Similar graphs can be drawn for other cases with available model parameters in Table 1.

The 'Low' discharges part in Figure 3 is represented by a straight line.

Other karstic aquifer and sinkhole discharge and hydraulic head measurements are available from the Newport (St. Marks River) region in the USA. Like the Crawfordville (Wakula River) case, its discharge and hydraulic head records are shown with three classifications as ‘High’, ‘Medium’, and ‘Low’. It is more obvious that in this case the recession limbs are extremely responsive with almost sudden falls. Again, the start and end coordinates of each partial straight-line limb segment coordinates are provided numerically in the second half of Table 1 with slope values, in addition to the model parameters c and d.

Applying Equations (7) and (8) to each discharge and hydraulic head parts of the records provides the graphs in Figures 5 and 6 for the discharge and hydraulic head values of the Newport karstic records.

Comparing these figures with Figures 2 and 3 indicates that in the Crawfordville (Wakula River) case there is a convex discharge recession, especially at the ‘V. high’, ‘High’, and ‘Medium’ sectors, but in Newport (St. Marks River) there is a concave reduction in the discharge recession part. In the latter case the recession is rather sharp without any significant fluctuations around the straight lines, which indicates that in Newport (St. Marks River) the aquifer is less karstifized than in Crawfordville (Wakula River). Furthermore, in Newport (St. Marks River), the sinkhole has four parts, each with hydrogeological features that act as a natural pipe sequence. This is due to insignificant contributions from adjacent fractures, fissures, and lateral solution cavities to the sinkholes. Finally, the comparison shows that the karstic region in Newport (St. Marks River) seems to be younger than the Crawfordville (Wakula River) aquifer.

## CONCLUSION

Karstic aquifers and sinkholes are spread all over the world, each with geographical and hydrogeological significant differences. They are non-homogeneous, anisotropic, and have a complex network of solution cavities, fractures, fissures and, at places, porous patches. In any theoretical analytical solution, one must depend on continuity and groundwater movement mathematical expressions to reach a general groundwater balance equation in differential equation forms. Although analytical approaches are valid mathematically, their adaptation to a given regional karstic aquifer of sinkhole mechanism is rather difficult and, therefore, for practical solutions further assumptions must be added. In this paper, rather than theoretical analytical holistic solutions, record pattern-based expert view solutions are presented based on partial verbal and then empirical mathematical solutions. The karstic domain discharge and hydraulic head records are subdivided into four parts as ‘Very high’, ‘High’, ‘Medium’, and ‘Low’ sections and each section is examined physically leading to practical solutions in empirical formulations. These formulations are valid universally, but their parameters must be determined based on the discharge or hydraulic head records in sinkholes if the medium is of karstic type. The application is presented for Crawfordville (Wakula River) and Newport (St. Marks River) karstic aquifer records, which are available in the literature. Instead of holistic mathematical solutions, piecewise solutions are presented with physical, logical, and practical views. It is recommended that in particularly karstic regions, sinkhole discharge and hydraulic head records availability must first be examined visually in a rational manner prior to discharge and hydraulic record holistic analytical modeling. It is recommended to have physical and practical interpretations with simple empirical formulations.

## CONFLICT OF INTEREST

None.