Modelling size constraints on carbonate platform formation in groundwater upwelling zones

Carbonate depositional systems related to groundwater upwelling are ubiquitous around the world and form ecologically and culturally important features of many landscapes. Spring carbonate deposits record past climatic and hydrological conditions. The reconstruction of past processes using spring carbonate proxies requires fundamental understanding of the factors that control their geometry. In this work, we show that the spatial extent of spring carbonate platforms is amenable to quantitative prediction by simulating the early growth stage of their formation for the iconic mound springs in the central Australian outback. We exploit their well-defined, circular geometry to demonstrate the existence of two size-limiting regimes: one controlled by the spring flow rate and the other by the concentration of lattice ions. Deviations between modelled and observed size metrics are attributable to diminishing spring flow rates since formation, enabling assessment of the relative vulnerability of springs to further hydrological change.

Stable oxygen and carbon isotope analysis in combination with dating techniques of fossiliferous carbonate deposits in groundwater upwelling zones have previously been used to derive changes in climate and the onset of arid periods, as well as information about the long-term temporal dynamics of groundwater upwelling [1][2][3][4][5][6][7] . With the exception of some works that have modelled the shape of actively-forming depositional features [8][9][10][11] , most previous studies concerning the geochemical and sedimentological controls on spring carbonate deposition are field-based, with a particular emphasis on understanding carbonate precipitation triggers 12 . The principal controls on the morphology of carbonate mounds have been assumed to be either the availability of lattice ions 13 or flow rates 14,15 . In this study, we used reactive transport modelling, i.e., numerical modelling of water flow and concomitant chemical reactions, to investigate the principal controls on the size of carbonate mounds associated with springs in the Kati Thanda-Lake Eyre South region of South Australia (Fig. 1). The springs in this region sustain unique wetland ecosystems and are of great cultural significance 15,16 ; consequently, methodologies that provide insight into their life cycle are needed for their management and to aid conservation efforts.
Mound formation is a three-part process consisting of: (i) initial spring and wetland formation, (ii) mound growth and (iii) growth cessation and stabilisation. Kinetically-controlled CO 2 degassing from emergent spring water as water flows away from the spring vent is central to this mound spring conceptual model. After a certain distance, calcium carbonate (CaCO 3 ) precipitation, triggered by CO 2 degassing, results in the build-up of a mound barrage that forms an enclosure for a central pool. The time required for sufficient degassing to trigger CaCO 3 precipitation determines both the size of the central pool and surrounding mound. The CO 2 degassing rate is dependent on flow turbulence and the surface area exposure of water to the atmosphere. Once a pool has formed, the degassing rate slows as the mound height increases and the water trapped in the central pool deepens and becomes less turbulent. Mounds reach maturity when the capacity of emergent spring water trapped within the spring pool changes from one of net precipitation to net dissolution and erosion, eventually leading to the development of a tail gutter through which water escapes to form a wetland delta at the mound base, to which carbonate precipitation shifts. This model not only explains many of the morphological and geological features of carbonate mound springs found in the Kati Thanda-Lake Eyre South region, but also the great ages of mounds, which may reach tens of thousands of years 16 relative to the comparatively fast precipitation rates found in the field 12 . Consequently, the early growth stage discussed in this paper describes a relatively short period of time in the life of a mound spring, but nevertheless an important one as it can provide a snapshot of the possible hydrological and hydrogeological conditions of a spring wetland and supporting aquifer at a given point in time.
Modelling past spring environments involves making assumptions about past hydrological conditions. Therefore, rather than calibrating a model to simulate observed mound sizes, this study's objective was to test if, given the most defensible set of assumptions of past flow conditions, a reactive transport model can provide: (i) an explanation for the observed large variation in mound radii, and (ii) quantitative support for our proposed conceptual model of mound formation 12 .
The studied pools, which are near-circular in shape, have a surface area between 100 and 420 m 2 , and associated well-vegetated wetlands. The mounds surrounding them are generally dome-shaped or shield-like (Fig. 2), with heights varying from 2 to 5 m. Their near-circular shape suggests that they initially formed by carbonate precipitation from discharging groundwater that flowed away from the spring vent in a radial fashion. This assumption is justified by the flat topography of the terrain. Consequently, the geometry of the earliest stage surface manifestation of a mound spring was conceptualised as a circular wetland (Fig. 3), with the discharging water radially flowing away from the vent.
Two end-member regimes are proposed to aid understanding of the factors limiting mound size. The first possibility is that the maximum radius of the mound is commensurate with the maximum distance that emergent spring water can flow away from the vent, which is limited by water loss in the surrounding wetland by evapotranspiration and infiltration. The second possibility is that the lattice ions (calcium and carbonate) that build the mound become depleted before the maximum migratory distance of water is reached, resulting in a mound radius smaller than the spring-supported wetland. For both regimes, a characteristic timescale exists. The timescale for the water-limited regime, t w , can be defined as the time required for the water to reach the edge of the wetland. The maximum extent of the circular wetland was based on water balance calculations, using measured spring discharges, as well as evapotranspiration and infiltration rates 13 (See Methods). As the calculated flow velocity becomes almost negligible towards the edge of the wetland, the t w was practically defined as the time require to reach 95% of the wetland radius. Changes in the water composition in the modern spring wetland in the study region have been successfully simulated previously by a reactive transport model that includes degassing of dissolved CO 2 and calcium carbonate precipitation 13 . With this model, a characteristic timescale for the lattice-ion limited regime can be determined. Here, this timescale, denoted t c , was taken as the time after which negligible calcium carbonate  . Graphical representation of the geometry and hydrological processes adopted for the modelling. R is the maximum radius of the spring wetland, r is the radial distance from the spring vent, Q is the spring discharge, ET and In are water-loss through evaporation and infiltration respectively, and h c is the mean height of the water column. precipitation occurs (see Methods). The carbonate precipitation rate is a complex function of the composition and pH of the spring water, which vary as degassing and carbonate precipitation reactions progress; the characteristic timescale must therefore be evaluated numerically. The analysis of timescales provides fundamental insight into the key factors that limit the extent of the spring mounds, but it does not provide any spatial or temporal constraints on the size of the mounds. Therefore, the current study also undertook reactive transport modelling to test if the size of the carbonate mound footprint, as well as the footprint of the central spring pool can be predicted based on knowledge of the spring discharge, water loss by evaporation and infiltration, and spring water composition. Consistent with the analysis of characteristic timescales, simulations considered CO 2 degassing and CaCO 3 precipitation, but within a spatially-explicit framework of radial-symmetric flow with water loss via evapotranspiration and infiltration. Such insights into the controls on carbonate mound morphology gained from the current investigation will also impact our understanding the palaeohydrology of the wetland and the supplying aquifer and therefore develop also our appreciation of the relative vulnerability of these wetlands and the associated groundwater resource to changes, both natural and anthropological.

Results
Spatial measurements taken from sampled springs in the Kati Thanda-Lake Eyre South region of South Australia are summarised in Table 1. The waters collected from the springs were brackish Na-Cl-SO 4 -HCO 3 -type waters. The temperature of spring waters varied between 19.4 °C and 30.4 °C and electrical conductivity (EC) ranged between 2.5 and 7.9 mS/cm. the Oxidation/Reduction Potential (ORP) indicated that with the exception of Blanche Cup and The Bubbler, spring waters were generally reducing in character upon emergence. All alkalinities were greater than 400 mg/L (as CaCO 3 ) and the pH of spring water was generally near-neutral on emergence (6.9 < pH < 7.4), apart from Beresford Spring (pH = 7.9). The partial pressure of CO 2 (P CO2 ) values for all emergent spring water were several orders of magnitude greater than average atmospheric values, and the saturation index for calcite (SI c ) ranged between −0.46 < SI c < 0.76. Although the water composition varies between springs, the Ca 2+ concentration is inversely correlated to alkalinity, and this relationship is controlled by the subsurface equilibrium of the water with calcite at a partial CO 2 gas pressure of P CO2 = 10 −1.5 atm (Fig. 4).
By virtue of this common control on the composition of the spring waters, a relationship between t c and the spring water Ca 2+ concentration exists that applies to the springs examined. By varying the Ca 2+ concentration (m Ca ) in the initial solution between 7.92 × 10 −4 <m Ca2 +<3.15 × 10 −3 mol/kg H 2 O, a relationship between the initial Ca 2+ concentration and the t c was defined by regression analyses, and the following dependence between t c and m Ca2+ was thus obtained (R² = 0.998): Equation (1) was used to draw the curved surface in Fig. 5, which marks where t w and t c are equal. Springs were classified as either water limited or lattice-ion limited based on their discharge and water loss rate (which determine t w ) and Ca 2+ concentration (which controls t c ). The data points for all but one of the five springs plot below this surface, indicating that lattice-ion availability is limiting mound size. Billa Kalina spring is the only spring that falls in the water-limited regime.
The simulated radius is compared to the mapped geological and morphological spring mound for each site in Fig. 6. Without calibration, the model adequately approximates the observed radius for three of the five springs (Fig. 7, Table 2). Another metric calculated by the model is the distance from the vent at which the maximum carbonate mineral precipitation rate occurs. At this point, the vertical growth of the mound is fastest, which ultimately results in the build-up of a rim behind which the central pool forms. Comparison of this distance with the size of the modern pool indeed shows some agreement of the model with the field data (Fig. 7). The vertical carbonate growth rates calculated by the model are in the order of millimetres per year, which is consistent with previously published experimental work 12 , as well as hydrochemical mass balance calculations 15 .

Discussion
Our findings significantly reveal that carbonate mound morphology is not only dependent on the spring's physical hydrological characteristics, but also, quantifiably, on the chemical composition of the emerging water. It explains why there is no correlation between these mound sizes and spring discharge; mounds of a high-discharge spring like The Bubbler remain relatively small as, due to the emerging water's chemical composition, fast precipitation kinetics ensures that most carbonate precipitation occurs close to the vent. The calculated vertical growth rates of a few mm/yr imply that the time required for mound formation is in the order of centuries to a few millennia 12 , which is short relative to the age of many mound springs 16 . Nearly all mounds in the Kati Thanda-Lake

Spring
Abr. Eyre South region are mature, thus it seems that the birth of new springs is infrequent at this period in geological time.
The model significantly underestimates the mapped radius of two of the four lattice-ion-limited mounds (Blanche Cup and Beresford), which are out by a factor of 3. This finding was expected though, as the present-day spring discharge used in the model is not necessarily representative for the discharge at the time of initial formation. In the Kati Thanda-Lake Eyre South region, spring discharge is believed to decrease with time 12,17 due to increasing mound height 14,15,17 , decreasing groundwater pressures caused by abstraction or long-term climate-related recharge rates 12 , and/or vent blockage 13 . By attributing the differences in the simulated and observed radii to a decrease in spring flow rates, the model can be used to provide an estimate of the spring  discharge at the initial stage of a mound's formation by adjusting the discharge to match the modelled and measured radius (Fig. 7). The resulting discharges (Table 2) remain well within the range of discharges of other springs in the region 18 . An increase of the present-day discharge by a factor of 25 is required to increase the simulated size of the Beresford spring mound to attain the observed size. The number of studied springs was necessarily limited in our study as only five sites had adequate aerial photographic coverage, minor disturbance from pastoralism and reliable recent discharge and hydrochemistry data in this remote region. Nonetheless, by virtue of the well-defined geometry of the spring mounds, our work clearly unveils the relationships between spring water chemistry, spring flow rate and wetland hydrology. Our findings equally apply to more complex carbonate deposits in which these fundamental controls are also expected to exist but are more difficult to recognise. For instance, similar-looking structures have been identified on Mars and their manifest morphological similarity has led to comparable hypotheses concerning formation via groundwater discharge 19,20 . Ancient martian hot springs have also been studied in Gusev crater by the Spirit rover mission 21 , strengthening the likelihood of spring mounds on that planet. The quantification of the size metrics of the mound structures by reactive transport modelling was consistent with field observations, while deviations are attributable to the known decline of the spring flow rates with time. This demonstrates the value of the methodology in conceptual model development and validation, in palaeo-environmental studies based on carbonate archives or in hydrological studies. Finally, its potential application, with appropriate modification, in remote areas and possibly even extra-terrestrial hydrology [19][20][21][22] where observations are largely limited to remotely-sensed data, potentially involving fluids other than water, is an exciting prospect.

Methods
Spring mounds at five locations within the remote Kati Thanda-Lake Eyre South region of South Australia were analysed as part of this study. The spatial extent of the spring mounds and the features of the surrounding terrain were mapped using stereoscopic pairs of aerial photographs. The exception to this was Billa Kalina Spring, where only a single image was available; however, the prominence of the mound compared to the surrounding landscape at this site made it clearly identifiable. All interpretations were verified in the field by tracing the perimeter with a GPS device. Surface area was determined using GIS software, and the effective radius of each mapped mound was calculated as if the surface area was circular.  A table of water quality measurements and chemical analyses collected for each of the spring localities in 2008 and 2009 is provided as Table 3. The field methods used to obtain water samples have been published previously 13 .
The characteristic timescale for water flow, t w , is defined as: where r is the radial distance from the spring vent (m), R is the maximum radial extent of the wetland (m), v w is the water flow velocity (m/s), which decreases with increasing r, and t w denotes the time it takes for the water to reach a distance of 0.95 R.
Steady-state conditions were assumed, i.e., the volumetric discharge rate of the spring Q (m 3 /s) equals the total rate of water lost via evapotranspiration plus the infiltration over the surface area A of the circular wetland. The value of v w as a function of r then follows from: where h c is the wetland water depth, e is the evapotranspiration flux (m/s), and i the water infiltration rate per unit of wetland area (m/s). The water-loss rates e and i were assumed to be spatially uniform and represent an average between warmer and cooler months.
In order to determine h c , analogous to modern spring environments in the area, it was assumed that after initial mound formation, the resultant wetland was colonised by dense sedge and reed vegetation quickly compared to the time estimated to build up the mound 15 . Flow within a shallow, vegetated wetland may be described as slow where q is discharge per unit width (m 2 /s), K w is the hydraulic conductance coefficient for overland flow (m 2−b /s), h r is the height of the surface water column at radius r (m), b is an exponent related to vegetation microtopography, stem depth and density distribution, S is the hydraulic gradient (dh r /dr), a is an exponent that expresses the degree of flow turbulence 26 . Values of 1 and 3 for the exponents a and b, and a K w value of 1.16 × 10 2 /m/s were adopted based upon previously published field-based experiments and subsequent recommendations for the modelling of wetlands with laminar flow and dense vegetation characteristics, such as those observed in mound spring wetlands 23,25,27 . Solving equation (4) for h r subject to the following boundary conditions: Equation (7) was numerically integrated with respect to r and the result was divided by R to find the average water depth (h c ) of the wetland that was used for subsequent calculations with Equation (3).
Discharge data to constrain Q were obtained from saline dilution tests or from temporary weir measurements over a series of field trips between October 2008 and July 2009. The evapotranspiration flux e and infiltration flux i were constrained using chloride and δ 18 O mass balance calculations 13 . Total water loss per unit surface area (e + i) was calculated by dividing the measured discharge by the area of the modern wetland environment, and was partitioned into rates for e and i using their relative proportions based on previously published data 13 . In this study, 33% of water loss was attributed to evapotranspiration and 67% to infiltration, which represents the average between results from warmer and cooler months.
Based on the reactive transport model by Keppel et al. 13 for modern spring environments, the changing chemical composition of the water as it flows away from the spring vent was conceptualised to be due to: (i) CO 2 degassing, (ii) evapotranspiration, and (iii) CaCO 3 precipitation or dissolution. PHREEQC-2 28 was used for modelling. This model was modified to reflect the flow conditions of the earliest stage of a spring, i.e., its initial surface manifestation.
The loss of dissolved CO 2 to the atmosphere by degassing was described using the following rate expression 29 :  13 for the model of the modern spring environment. This value is also within the range of published values for stream environments [30][31][32] . Carbonate precipitation was based on the rate expression for calcite dissolution by using the kinetic rate equation by Plummer et al. 33 , which is implemented in the WATEQ. 4 F database 34 of PHREEQC-2.
The modelling to quantify the timescale t c was conducted using PHREEQC-2 by simulating CO 2 degassing and kinetically-controlled carbonate precipitation in a batch-type model. Initial solutions were defined by equilibrating the solution with calcite and CO 2 gas at a partial pressure of 10 −1.1 atm. The latter value was chosen based on the relationship between alkalinity and the concentration of Ca 2+ according to the field data (Fig. 4). Sodium (4.688 × 10 −2 < m Na < 6.416 × 10 −2 mol/kg H 2 O, concentration m Na was varied to achieve charge balance) and Cl − (5.093 × 10 −2 mol/kg H 2 O) were included as background electrolytes to obtain the approximate ionic strength of the spring waters. Calcite precipitation was considered to be negligible when the calculated precipitation rate  Table 3. Water quality and chemistry results for the spring sites included in this study. All samples were taken from the spring vent. a All analytical determinations in mg/L unless otherwise stated.
SCieNTifiC REPORTS | (2018) 8:17460 | DOI:10.1038/s41598-018-35771-z dropped below 10 −10 mol/kg H 2 O/s, which is roughly 3 orders of magnitude lower than the maximum rate found in the simulations. PHREEQC-2 uses a mixing-cell approach to simulate transport along a flow-line, and the flow velocity is defined implicitly based on a fixed transport time step length and the length of a grid cell. Cell lengths were calculated by substituting v r = dr/dt into Equation (3)  In which n denotes the cell number, α π = + e i ( ) and β = 2πh c . Solving for r n+1 , gives: By setting a time step (t n −t n+1 ) and assigning a radial distance (r 1 ) for the first cell, equation (10) was used to calculate the lengths of the consecutive cells. The one-dimensional discretisation thus varied for each of the five springs modelled, depending on their flow rate and the evapotranspiration and infiltration fluxes ( Table 1). The size of the mound was taken as the distance along the flow path where the calcite precipitation rate fell below 10 −10 mol/kg H 2 O/s, which is less than 0.5% of the maximum rate observed (Fig. 8).