Estimation of deep infiltration in unsaturated limestone environments using cave lidar and drip count data

Limestone aeolianites constitute karstic aquifers covering much of the western and southern Australian coastal fringe. They are a key groundwater resource for a range of industries such as winery and tourism, and provide important ecosystem services such as habitat for stygofauna. Moreover, recharge estimation is important for understanding the water cycle, for contaminant transport, for water management, and for stalagmite-based paleoclimate reconstructions. Caves offer a natural inception point to observe both the long-term groundwater recharge and the preferential movement of water through the unsaturated zone of such limestone. With the availability of automated drip rate logging systems and remote sensing techniques, it is now possible to deploy the combination of these methods for largerscale studies of infiltration processes within a cave. In this study, we utilize a spatial survey of automated cave drip monitoring in two large chambers of Golgotha Cave, southwestern Western Australia (SWWA), with the aim of better understanding infiltration water movement and the relationship between infiltration, stalactite morphology, and unsaturated zone recharge. By applying morphological analysis of ceiling features from Terrestrial LiDAR (T-LiDAR) data, coupled with drip time series and climate data from 2012 to 2014, we demonstrate the nature of the relationships between infiltration through fractures in the limestone and unsaturated zone recharge. Similarities between drip rate time series are interpreted in terms of flow patterns, cave chamber morphology, and lithology. Moreover, we develop a new technique to estimate recharge in large-scale caves, engaging flow classification to determine the cave ceiling area covered by each flow category and drip data for the entire observation period, to calculate the total volume of cave discharge. This new technique can be applied to other cave sites to identify highly focussed areas of recharge and can help to better estimate the total recharge volume.


Introduction
Karstic aquifers represent substantial global groundwater resources (Worthington and Gunn, 2009).However, the phenomena related to water movement in the unsaturated zone of karstic systems are not yet fully understood.To better manage karst resources, it is important to understand and predict how water flows in karstified limestone.Many traditional methods developed for modelling groundwater flow regimes in highly heterogeneous karstic aquifers are focussed on the faster drainage components, i.e. conduits and channels (Morales et al., 2010(Morales et al., , 2007;;Pardo-Iguzquiza et al., 2011;Smith et al., 2012;Ford and Williams, 2007;Goldscheider and Drew, 2007).However, these methods are less suitable in characterizing water movement through the smaller fracture or matrix flow components of the unsaturated zone, lacking vital information relevant to the complete understanding of flow through fractured rocks.The formation of combined flow networks is the key phenomenon that separates karst aquifers from porous and fractured-rock aquifers (Ghasemizadeh et al., 2012).
Recharge estimation is critical for understanding the water cycle and contaminant transport, and for water management.However, monitoring water in the unsaturated zone, especially in highly heterogeneous limestone formations, is difficult, and estimating recharge is one of the most complicated tasks in the hydrological cycle (Scanlon, 2013).Continuous water content measurement using time domain reflectometry (TDR) (Rimon et al., 2007;Dahan et al., 2007) or neutron activation (Koons and Helmke, 1978;Sophocleous, 1991) allow point study on the unsaturated zone water infiltration rate.Tracers such as fluorescent dyes and environmental isotopes in the unsaturated zone at many sites showed an order of magnitude range in recharge rates over 7-70 m yr −1 (Sheffer et al., 2011).This has been attributed to different flow systems (quick flow and slow flow), arid versus humid climate forcing, and variations in storage in the soil and epikarst, e.g.Mendip Hills, England (Friederich and Smart, 1982), Israel (Even et al., 1986), Niaux, France (Bakalowicz and Jusserand, 1987), Pennine karst, England (Bottrell and Atkinson, 1992), Slovenia (Kogovšek, 1997), and Mt Carmel, Israel (Arbel et al., 2008).
At the scale of cave drip waters, studies in geologically old, fractured limestone (that has undergone past diagenesis) have identified the importance of fracture flow and storage in solutionally enhanced fractures or caves.Matrix storage is also possible in geologically young limestone that has not undergone past diagenesis and contains primary porosity.Cave drip waters are fed directly from the karst bedrock and overlying soil (White, 2002;Tooth and Fairchild, 2003;Atkinson, 1977;Raesi and Karami, 1997;Ford and Williams, 2007), and variations in the size and orientation of fractures, together with variable storage capacity, play fundamental roles in governing the response of a drip site to individual recharge events (Baker and Brunsdon, 2003).High secondary porosity is associated with the epikarst, a zone of heavily weathered carbonate rock, which may act as a water storage reservoir retarding flow and sustaining slow percolation through the unsaturated zone in rocks where karstification has occurred (Arbel et al., 2010;Williams, 1983).
Drip discharges have been categorized in terms of the type of flow process occurring between recharge water and drip water.One possibility is that the drip water is transported via direct delivery of recharge along preferential flow paths (e.g.fractures).Another one is piston flow, where stored water is expelled from pores and fissures by incoming infiltration water (Baker et al., 2000;Tooth and Fairchild, 2003).A more refined understanding of karst infiltration has been achieved through the use of continuously recording (automated) drip measurement devices, which are capable of resolving fine temporal changes in drip rate (McDonald and Drysdale, 2007).Studies incorporating such measurements have increased our knowledge of seepage dynamics.For example, Markowska et al. (2015) classified five different drip types at Harrie Wood Cave in south-eastern Australia, suggesting the heterogeneous flow in the unsaturated zone due to the nature of the karst architecture.Jex et al. (2012) classified drip behaviour using high-resolution drip time series and employed multi-dimensional scaling (MDS) analyses to cluster the data accordingly.Studies using automated counters have also discovered the role of atmospheric pressure on drip variation, and questioned the linearity of rechargedischarge response at various timescales (Genty and Deflandre, 1998;Baker and Brunsdon, 2003).
Caves offer a natural inception point to observe both the long-term recharge and the preferential movement of water through the unsaturated zone of such fractured bedrock by monitoring stalactite drip rates.With the availability of both new drip rate logging systems and remote sensing techniques, it is now possible to deploy the combination of these new measurement methods for larger-scale studies of many individual drips within a cave.The goal of this paper is to demonstrate the nature of the relationship between flow types classified by the morphological analysis of stalactites (Mahmud et al., 2015) and drip time series characteristics.A spatial survey of automated cave drip monitoring in two large chambers of Golgotha Cave, south-western Western Australia (SWWA), is utilized to achieve this goal.Recharge into the cave is quantified based on the drip data, Terrestrial Li-DAR (T-LiDAR) measurements and flow classification.Li-DAR is a word which combines "light" and "radar", although the word LiDAR is thought by some to be an acronym for light detection and ranging.We estimate the water balance to develop a simple model describing the ground surface extent from which flow is focussed on the monitored cave ceiling area and the associated lateral flow within the Tamala limestone formation.

Site description
Our study site, Golgotha Cave (36.10 • S, 115.05 • E, Fig. 1a), is developed in the Spearwood System of the Tamala limestone, which comprises medium to coarse-grained Quaternary calcarenites of predominantly aeolian origin.The limestone is wind-blown calcareous sands that have deposited widely around the western and southern coasts of Australia.The cave chamber was formed by unsaturated zone water flow and subsequent widening by ceiling collapse.The cave is 200 m long and up to 25 m wide, and the limestone bed is 20-30 m thick over the cave (Fig. 1b).The surface vegetation is wet eucalypt forest with a substrate of weathered siliceous dune sands.The porous sand layer overlies the dune limestone, to variable depths, measured up to 3 m deep (Treble et al., 2013), with more than 30 m underlying dune limestone.The overlying surface of the Tamala limestone is mantled by sands.Therefore, few karst features are seen on the ground surface except for the occasional outcrop, and the presence of these are far fewer than the conventional karst landscape.There is no karren within this karst landform.There are occasional dolines, e.g. the cave entrance, formed by cavern col-  Mahmud et al., 2015).Chamber 1 contains sampling area 1 (T-LiDAR 1) and chamber 2 contains sampling area 2 (T-LiDAR 2) of the drip water chemistry monitoring programme operating since 2005 (Treble et al., 2015(Treble et al., , 2013)).
lapse, but none above the study sites.Additionally, in chamber 2 there is evidence of solutional widening generally associated with high-flow sites (Treble et al., 2013).The overall land surface gradient above the cave site is calculated as ∼ 20 % using the topography map which has the cave outline properly registered on it.However, there is negligible surface runoff given the high permeability of the sand layer, with hydraulic conductivity ranging from 100 to 2000 m day −1 (Smith et al., 2012).Over several field campaigns, runoff has never been observed at the site.This indicates zero surface runoff and maximum potential for infiltration through this karstified limestone.
SWWA has a Mediterranean-type climate, with warm to hot, dry summers (mid-November to mid-February) and mild to cool, wet winters (mid-May to mid-August), associated with the seasonal migration of the mid-latitude westerly winds.Rainfall recorded since 1926 at Forest Grove is 1136.8± 184 mm annually (BoM, 2015) (34.07 • S, 115.10 • E, weather station number 9547, 5 km from the site; Fig. 1a).Typically, the highest rainfall starts in late autumn (May) and carries on during the entire Southern Hemisphere winter wet period (May-October) (median monthly rainfall is approx.100 mm) (Fig. 2a).Mean maximum daily temperatures range from 16 • C in July to 27 • C in February.Moisture is delivered by troughs embedded in the westerly flow, but may be sourced from regions to the south-west or northwest of our site (Bates et al., 2008;Fischer and Treble, 2008).Recorded monthly rainfall conditions are shown in Fig. 2a, for hydrological years 2012, 2013, and 2014, starting from April when the water budget is close to zero.Each hydrological year has a similar pattern during the dry period, with months from October to April showing a water deficit or only a negligible amount of recharge (Fig. 2b).In contrast, there is significant amount of excess water available to infiltrate during the wet season.The hydrological year 2014 was rather dry, having 943.8 mm precipitation, far below the long-term annual mean precipitation (1141 mm).Year 2013 was relatively wet (1239.8mm), whereas 2012 (1008.6 mm) was slightly below the long-term annual mean.
Cumulative water budgets are calculated using precipitation (P ) and modelled evapotranspiration (ET) data from the Australian Water Availability Project (AWAP) (Raupach et al., 2009) in order to obtain the total infiltration in the karstic aquifer for each hydrological year.According to Raupach et al. (2009), the ET is modelled using the mathematical equations of Priestley-Taylor (Smith et al., 2012).Monthly calculated evapotranspiration is subtracted from the monthly rainfall totals to determine the water budget (i.e.P −ET) shown in Fig. 2b.Potential infiltration is then calculated from all positive monthly water budgets (monthly excess water).The total sum of all monthly excess water for a hydrological year (from April to March the following year) gives the potential infiltration for that particular year.We calculate an average annual potential infiltration of 696.3 mm during 3 observed hydrological years, ranging from a minimum of 608.6 mm in 2014 to a maximum of 858.7 mm in 2013.
We installed 34 drip water monitoring sites in the wettest areas of two large chambers of the cave (Fig. 1b) named Li-DAR 1 (i.e.sampling area 1, located approx.60 m into the cave) and LiDAR 2 (i.e.sampling area 2, located approx.90 m into the cave).Figure 3 shows the studied ceiling area above the loggers in each chamber.The notation used for site identification consists of an Arabic numeral and a letter/Roman numeral.The first Arabic numeral indicates the chamber and the following letter/Roman numeral indicates a certain drip site within the given chamber.In the second position, a letter/Roman numeral is assigned to distinguish between drip data collection processes: the letters specify the sites having both manual and automatic drip counts, and the sites with Roman numerals only have drip logger data.In each chamber, the drip loggers were laid out in rough transects approximately following the ceiling gradient.In chamber 1, sites 1A, 1B, 1i, 1ii, 1iii, 1iv, and 1v are underneath the northern side of the ceiling with slightly thinner limestone overburden (32 m thick), compared to the southern side (32.6 m overburden) where drip loggers 1vi to 1xi were placed (Fig. 3a).This variation in overburden thickness within chamber 1 represents an overall 0.6 m lowering in ceiling elevation from the northern side to the southern side.This slight variation in ceiling elevation means that a higher hydraulic gradient occurs at the southern patch that is more densely decorated with stalactites.In chamber 2, the sites are spread over a larger area.The southern portion is close to the intersection of the ceiling with the cave wall, having comparatively low ceiling elevation and high overburden thickness.On the other hand, the northern side of Fig. 3b is far away from the wall.Site 2E is located in the wettest area, close to the lowest point at which the ceiling and wall intersect, whilst 2B is ∼ 5 m from the wall and 2A is ∼ 10 m from the wall.

Background on flow type classification
The literature suggests that karst hydrological flow properties can be identified from the geometry of stalactites and other morphological features in relation to the cave ceiling (Fairchild and Baker, 2012).Based on this concept, Mahmud et al. (2015) used T-LiDAR measurement to image a cave ceiling including individual stalactites.Statistical and morphological analyses of the point clouds produced by T-LiDAR were then used to categorize the ceiling features into different flow types.Through this methodology, the role of the type of water flow processes was analysed and a preliminary conceptual model was developed by studying the spatial distribution of a large population of stalactites and their geometric properties (length and aspect ratio) in three sites within the Golgotha Cave system.

Flow type classification based on LiDAR data
In this section we briefly describe the methodology that was used to investigate different flow patterns classified in Mahmud et al. (2015).Based on the typical types of porosity and infiltration processes in karst, Mahmud et al. (2015) defined three categories of flow for the observed ceiling morphological features.These are matrix flow, fracture flow, and a combination of conduit, fracture, and matrix flow.The matrix flow category was further subdivided into two subclasses: soda-straw and icicle-shaped stalactites based on their geometric properties.The LiDAR-based morphological analysis was used to identify individual stalactites and flow classifications (icicle, soda straw, fracture, and combined flow) based on the T-LiDAR point clouds.Figure 4a shows the topography (in 2-D) of a portion of the ceiling at site 1 classified into the different flow categories.The outcome of this classifica-tion is shown in Fig. 4b. Figure 5 shows a similar analysis for a different ceiling portion of site 1 containing stalactites feeding sites 1ii, 1iii and 1v classified as soda straw (Fig. 4c).
According to the morphological analysis shown in Figs. 4 and 5, stalactites feeding sites 1A, 1B, 1i, 1ii, and 1iii are in the icicle flow category, and the stalactite feeding site 1v is classified as soda straw.Similar morphological analyses were performed for all stalactites that feed the loggers in both studied chambers.As a result, we identify that the majority of drip flows from chamber 1 are icicle type, with one possible location for each of soda straw (site 1v), fracture (site 1x), and combined flow (site 1viii).In contrast, drip loggers in chamber 2 record a completely different setting, having a variety of flow patterns (nine locations for icicle flow, four for fracture flow, three for combined flow, and two soda straws).Based on the preliminary conceptual model developed by Mahmud et al. (2015), it was determined that chamber 1 is controlled by a matrix flow pattern characterized by icicleshape and soda-straw stalactites that are widely distributed in the roof of chamber 1 (Fig. 3a).Within such a system, infiltration rates are directly proportional to the matrix permeability representing the primary porosity of the karst formation.Hence rates of change of water movement are expected to be low, with slow drip rates of low variability (Fairchild and Baker, 2012).The stalactite pattern in the ceiling of chamber 2 is shown in Fig. 3b.Morphological analysis of Mahmud et al. (2015) showed that this chamber is likely to be controlled by fracture and combined flow; hence, drip rates are expected to vary over time in this chamber, depending on water transport in the preferential flow system.These results are consistent with field evidence: the chamber 1 ceiling is dominated by soda-straw and icicle-shape stalactites suggesting matrix flow, and stalactites in chamber 2 tend to be more isolated and evolving either along fractures or at the margins of relict dune surfaces (Treble et al., 2013).

Data acquisition and methodology
We investigate the relationship between infiltration through the fractured limestone and cave drip water discharge using the morphological analysis of ceiling features, coupled with drip logger and climate data from 2012 to 2014.In this paper, we locate the individual stalactites feeding the drip loggers using T-LiDAR images and digital photos, and identify each as matrix (soda straw or icicle), fracture, or combined-flow.These morphology-based classifications are compared with flow characteristics from the drip logger time series.The discharge from each stalactite is calculated and used to estimate the total discharge over each studied area.The total discharge from each area is compared with infiltration estimates to better understand flow from the surface to the cave ceiling of the studied area.

T-LiDAR and elevation data
Ground-based T-LiDAR is a commonly used remote sensing technology that records high-resolution three-dimensional point clouds of the Earth's surface, and its use in geology has been growing in recent years (Pringle et al., 2006;Fabuel-Perez et al., 2009;Rotevatn et al., 2009;Wilson et al., 2009).However, karstic model development using T-LiDAR is a novel application.There are few studies that discuss the benefits and the use of this tool, as well as the methods needed to work with this kind of technology (exceptions are Zlot andBosse, 2014b, a, andKaul et al., 2016).We have used a FARO Focus3D Terrestrial LiDAR to acquire threedimensional geological images of the cave ceiling, which is capable of capturing millions of three-dimensional points coordinates within a few minutes.The reason we used LiDAR was that it allows us to determine the precise position of each ceiling stalactite and its size and shape, which is very difficult to obtain with photographs or with any other remote sensing technique.Moreover, to locate the ground coordinates of individual drip loggers in the cave is also challenging, and a high accuracy is necessary to pair each logger with a stalactite.The advent of relatively low-cost LiDAR makes this study possible.The T-LiDAR measurements were taken adjacent to two locations where cave drip waters have been sampled for drip rate and chemistry for the past 10 years (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) (Treble et al., 2015(Treble et al., , 2013)).T-LiDAR positions 1 and 2 (Fig. 1) were selected such as to cover the significant portions of the ceiling from a perspective close to vertical in order to capture all 34 stalactites feeding the drip loggers, while minimizing the occlusion of other stalactites further away from the scanner line-of-sight.Figure 3 shows scan images for both chambers.
The elevation of the cave floor at the gate is 72 m a.s.l.according to GPS measurements.Elevations of the drip loggers inside the cave were obtained by a cave survey using a fibre surveyors tape as well as a SUUNTO tandem 360PC/360R clinometer to calculate dip and orientation (estimated error for dip is ±1-2 • , and for orientation it is ±5 • ).A metal surveyors tape was used to measure distance (accuracy ±0.1 m).Ceiling heights were surveyed from the T-LiDAR data that consist of both the drip logger sites and the ceiling features (i.e.stalactites).A detailed surface survey was performed to gain surface elevation and an estimate of the total thickness of overburden over the cave.A Bosch GLR225 Laser Distance Measurer was also used to measure the distance between points (accuracy ±1.5 mm).

Alignment of drip loggers and stalactites
To build the relationship between the flow patterns classified in Mahmud et al. (2015) and the drip time series characteristics, we initially start pairing the individual stalactites feeding the drip with the loggers time series.Close-up images from the ceiling areas were analysed to identify the exact stalactite locations by comparing T-LiDAR images, digital photos, and on-site observations.Figure 3 shows the studied ceiling area of both chambers captured using a T-LiDAR, with stalactite locations identified in blue circles.The stalactites (in short, Stal) are named according to the drip logger site (e.g.Stal 1A feeds the logger site 1A, and so on).Points (1, 2, 3, 4, 5) are assigned to some significant ceiling features (for example larger stalactites) to register the locations in the T-LiDAR images, digital photos, and on-site within the cave (Fig. 3).

Drip logger data
Automatic drip monitoring sites were established in August 2012.Stalagmate drip loggers (www.driptych.com)were installed throughout the two large chambers (Fig. 1).Data loggers were set to record the number of drops that fell on the logger per 15 min.Data were downloaded at regular intervals of 6 months, between August 2012 and March 2015, and data collection is ongoing.Preliminary screening of all drip time series was performed for quality assurance.Based on the initial data screening, we entirely discard five drip sites, i.e. 1iv, 1vii, 1xii, 2ii, and 2xii.Drips 1iv, 2ii, and 2xii show sudden changes in drip rate that likely reflect the logger being accidentally moved or misaligned after data downloads.Loggers 1vii and 1xii were discarded due to the recording of dual drips or missing data.The remaining 29 sites are con-sidered in this study, although parts of these time series were discarded where we considered the data unreliable.Data recorded during periods of known fieldwork were removed from the drip rate time series, including 1 day either side of recorded field trip days as standard protocol.Finally, the time series gaps are filled with synthetic data using spline interpolation of the drip data considering the drip statistics and correlation between drip rates.
The resulting drip rate time series are plotted for both cave chambers and 3 hydrological years from August 2012 to March 2015 in Fig. 6.Chamber 1 drip loggers show two different groupings in terms of flow rates: one group has drip rates of less than five drips per 15 min and another one has a higher drip rate.However, most chamber 1 drip loggers exhibit a clear response to the 2013 wet winter, presenting peaks at the end of September 2013 (Fig. 6a).On the con-trary, chamber 2 drip rates are more variable between sites (Fig. 6b-d).To clearly visualize the drip time series of chamber 2, we illustrate all these time series in three different plots based on their flow rates throughout the 3-year study period: (i) static drips with little discernible variation and very low flow rates (Fig. 6b), (ii) medium-variability drips with moderate discharges (Fig. 6c), and (iii) high-variability drips with high discharges (Fig. 6d).Each figure is scaled differently due to the differences in the magnitude of drip values.For easy identification of drip time series, we plot a few representative time series rather than all drips in Fig. 6b.Comparing the slow dripping sites from both monitored chambers, we discover a persistent base flow component even during periods of water deficit, which feeds storage water to the drip site.This indicates that the Tamala limestone formation does have significant storage within its matrix porosity to deliver uninterrupted cave discharge during dry periods.These monitoring sites (Fig. 6) exhibit a near-constant drip rate, with little or no relationship with hydrologically effective precipitation, except for the response to the 2013 wet winter.Nonetheless, there is still considerable variation between these static drips in terms of base flow, magnitude of response, and attenuation.Records presented in Fig. 6a (particularly site 1viii) and Fig. 6b (most of the slow drip time series) have noisy data with high variation over a short period.When drip rates are of approximately the same frequency as the logging interval (15 min), drip variability increases and are an artefact of the sampling interval.This is observed at site 1viii and most of the slow drip sites from chamber 2. The observed variability is typical of that of previous studies using drip loggers.We presented the raw data that are not corrected for barometric pressure variations -barometric effects explain the drip variability over the time period of days to weeks, and have previously been observed in porous limestone (Genty and Deflandre, 1998).

Cave discharge estimation
The infiltration through the limestone formation within the monitored areas of both chambers is estimated based on the drip data and the T-LiDAR measurements.We consider the drip data for the entire observation period to calculate the total drip counts for each logger and thus obtain the total cave-integrated drip water volume, considering 1 drip = 0.1433 mL according to Genty and Deflandre (1998).We identify 29 individual stalactites that relate to the logger drip data we placed in two large chambers of this cave site.We then extrapolate these logger data to the entirety of both chambers to predict the total infiltration within these areas.This is a fundamental outcome of this study which is made possible by the combination of drip loggers and LiDAR technology.The methodology does not include the area covered by each flow class; rather, it counts the total number of stalactites that fall into various flow categories.For example, in chamber 1 we find a total of 1909 individual stalactites in the ceiling; among these 1649 are considered to have matrix flow, 17 are fracture flow, 3 fall into a combined flow category and the remaining 240 are soda straws.Furthermore, we know the average drip discharge for each flow category (q), thus potentially allowing us to estimate the total flux volume.
Not all stalactites are actively dripping; therefore, to have an accurate quantification of the total drip flux, we need to know the fraction of stalactites that are actively dripping.It is not possible to identify the water beads hanging from active stalactites based on the LiDAR images.For this reason we cannot automatically count the active drips on large ceiling areas and we have to use photographs taken with flash that are only representative samples for each chamber, and we use an average value for each chamber individually, rather than for each flow category.Thanks to the flash, each water bead is visible as a shiny spot.A series of the digital photos of the chamber ceiling is analysed and the stalactites with a water bead on their tip are manually counted.One single frame from the chamber 1 ceiling is shown in Fig. 7.In this ceiling portion, there is a total of 45 stalactites (total green and yellow circles), among which 32 are actively dripping (green circles).Therefore, in this image 71.1 % of the stalactites are active.Similarly, we use other digital images of ceiling areas covering the rest of the monitored sites, and consider the average percentage for each chamber ceiling.For example, we apply the average proportion of actively dripping stalactites (57 %) to the total counts of different stalactites to only count the active stalactites.However, LiDAR just provides a onceoff snapshot and does not capture the temporal variability of actively dripping stalactites.Therefore, repeat surveys based on LiDAR are feasible (Zlot and Bosse, 2014b, a;Kaul et al., 2016) and would be a rapid and easy way of looking at changes over time.

Relation between LiDAR-based classification and drip data
The LiDAR-based flow classification for individual sites is detailed in Tables 1 and 2. The average drip rates per 15 min, skewness, coefficient of variation (COV), and the elevations of cave ceiling at the location of the stalactites are also listed to compare various flow categories.Coefficient of variation (COV) is defined as the ratio of the standard deviation (σ ) to the mean (µ).To build up the relationship between LiDARbased flow classification and drip rates, drip logger mean discharge, skewness, and COV are plotted in Fig. 8 against the overburden limestone thickness above all drip sites (from 32.4 to 41.9 m).We observe no significant relationship between the drip logger mean discharge, skewness, and COV with overburden thickness; however, these properties can be used to characterize different flow classes, and are discussed below in further detail.
All chamber 1 icicle flows have slow drips, i.e. less than five drips every 15 min (Table 1), the resulting mean discharges ranging from 6 to 20 L yr −1 .Their drip rates remain almost constant throughout the study period (Fig. 6a).However, it is evident that the drip loggers exhibit a clear response to the 2013 wet winter, with peak discharge at the end of September 2013 (Fig. 6a).Icicle drip discharges in chamber 2 show similar results, with slightly higher drip rates up to 8.5 drips per 15 min, discharging 42.7 L yr −1 .All these icicle flow types display low values of skewness and intermediary variation (COV) within the time series (Fig. 8b and  c).
Drip sites classified as soda straws (site 1v, 2iii, and 2xi) usually have very low discharges (less than two drips per 15 min).However, these drips have large skewness and variation (COV) (Tables 1, 2, and Fig. 8b).Among these soda-straw flows, site 2xi has the lowest mean discharge of 0.5 L yr −1 and shows the largest variation (COV) within the time series (Table 2).Moreover, such soda-straw drips do not exhibit constant discharge according to the drip time series (Fig. 6a and b), even though they have extremely low discharges.
Sites classified as combined flow type (1viii, 2E, 2v, and 2viii) have high discharges ranging from 12 to 28 drips per 15 min (60-140 L yr −1 ).In addition, these drips have a comparatively extended range of skewness and COV (Fig. 8b  and c).Lastly, sites 1x, 2i, 2vi, 2ix and 2xvi, characterized as fracture flows according to the morphological analysis, typically have the largest discharges (Fig. 8a).There is significant variability in discharge between these sites and within the individual time series (Fig. 6), evidenced by differences in mean discharge rates, skewness, and coefficient of variation, e.g.mean discharge ranging from 90 to 2700 L yr −1 .We observe rates of water movement ranging from 0.5 to 6.5 L yr −1 for soda-straw stalactites and up to 43 L yr −1 for icicle flow category, 60 to 140 L yr −1 for combined flow systems and 90 to 2700 L yr −1 for fracture flow type.This finding is different to recent studies in Mt Carmel Cave in Israel, with higher cave discharges of 320 mm h −1 ≈ 2.8 × 10 6 L yr −1 (Sheffer et al., 2011) and 1.9 to 3.5 × 10 6 L yr −1 (Arbel et al., 2010) for slow drip sites.The combined and fracture flow drip rates are also significantly lower compared to Mt Carmel Cave (2.8 × 10 6 L yr −1 for the intermediate flow to > 10.5 × 10 6 L yr −1 for quick flow (Arbel et al., 2010;Sheffer et al., 2011)).However, our drip discharge variations agree with other studies across Australian cave sites (Cuthbert et al., 2014;Markowska et al., 2015;Treble et al., 2013;Jex et al., 2012).
We compute the correlation matrix between all drip sites (Fig. 9).The different sites are aligned in the matrix according to the flow classification.Various flow types (I: icicle; F: fracture; C: combined; S: soda straw) characterized by the morphological analysis are shown in parentheses with the site notation in Fig. 9.The correlation between similar flow types can be visualized from the correlation matrix (Fig. 9).All chamber 1 sites from the northern patch representing icicle and soda-straw flows (1v, 1A, 1B, 1i, 1ii, and 1iii) are highly positively correlated (Fig. 9a).However, sites 1vi, 1ix, and 1xi (with icicle flow characteristics and falling within the southern patch of chamber 1) do not indicate significant correlations, suggesting spatial dependence on the drip discharge.On the other hand, chamber 2 sites belonging to the fracture and combined flow categories show moderate to high positive correlation with each other, possibly being highly responsive to rainfall events (Fig. 9b).All icicle and sodastraw flow sites are correlated with each other in a similar fashion, except sites 2A and 2xiii, which are negatively correlated with the rest of the drip discharges.These two sites show decreasing drip rates, while the rest of the sites from chamber 2 have increasing trends, suggesting flow-switching possibility at higher flows/heads.

Cave discharge
The total flux (Table 3) for each flow category i is obtained (i.e.Q icicle , Q fracture , Q combined , Q soda-straw ) by multiplying the total number of active stalactite for each category (n) by the corresponding mean measured drip discharge for the corresponding category (q):  Q combined = n combined q combined , Q soda-straw = n soda-straw q soda-straw . (1) The total discharge (Q) is obtained by summing the flux for all categories: We predict total cave discharge amount within the monitored ceiling area of both chambers using total drip counts and flow types categorized by the morphological analysis.Table 3 demonstrates the summary for both chambers considering all drip discharges.Based on our LiDAR and cave drip data, we have established a conceptual model for the cave site (Fig. 10).Figure 10a shows the vertical section of the entire cave with possible dune bedding and groundwater flows.In this study, we have estimated the water balance to understand the extent to which flow is focussed on the cave ceiling and to quantify the associated lateral flow within the Tamala limestone formation (Fig. 10).Monitored ceiling areas are 30.4 and 55.2 m 2 respectively for chambers 1 and 2, though the infiltration water comes from a larger surface area of 36.8 and 92.9 m 2 (Fig. 10b and c), suggesting that infiltration is being focussed on the studied areas of each chamber, by approximately 120 and 170 % respectively.In the conceptual model, we show near-vertical water movement, but there could also be lateral movement along the dune bedding.The bedding shown in Fig. 10 is based on the geometry that we observe in the lowest part of the cave, which is the third facies.For the upper part of the cave, the depicted bedding direction/angles are indicative only because it is much harder to see in the ceilings of chambers 1 and 2, due to the weathered surfaces and the ceilings in the studied areas that are relatively featureless.Therefore, the red column represented in Fig. 10b and c could be moved to fit the geomorphology and could change angle due to the dune bedding.However, a major finding from this research is that the cave seems to have a capillary fringe effect, with very little recharge entering the cave compared to the overlying surface infiltration.So lateral flow may be important, but it will likely also be affected by the capillary effect.A portion of this calculation may be the result of uncertainties involved in the methodology such as modelled ET or active stalactite count, but our field observations of the cave ceiling support focussed flow.For example, these studied ceiling areas are dominated by karst water infiltration compared to other ceiling portions consisting of dry and bare rock surface without any stalactite formation.Moreover, focussed diffuse flow is evidenced in Golgotha Cave by saturated rock viewed in vertical cross section in the cliff face and from within the cave ceiling, by clustering of sodastraw stalactites (Treble et al., 2013).Another possible uncertainty source involves the process of stalactite identification and flow classification based on the morphological analysis, which controls the amount of measured flow.This is a semiautomated segmentation process where the user controls parameters such as moving average window size and threshold to perform the morphology-based analysis.These parameters are adjusted to the particular cave site and, typically for segmentation procedures, involve a trade-off between over-and under-segmentation (Hyyppa et al., 2001).While the optimal parameters minimize both types of error, there is generally no parameter combination that results in a perfect, error-free segmentation (Mahmud et al., 2015).
The total counts of stalactites identified by morphological analysis in both chambers are around 2000, with chamber 1 being almost twice as densely populated as chamber 2 (Table 3).However, the proportion of active stalactites is higher in chamber 2, with a larger flux (per m 2 ) that potentially suggests greater lateral flow dispersion in chamber 2 (Table 3).Moreover, the chamber 2 ceiling covers a greater surface area (92.9 m 2 ), has the relatively lower ceiling elevation, and is adjacent to the cave wall (Fig. 10).This points to an influence of hydraulic gradient on water movement in this area.We suggest that this may indicate that a large portion of infiltrating water is flowing around the cave and inside the ceiling, rather than through it (Fig. 10c).The cave ceiling may be acting as a capillary barrier resulting in water moving along the ceiling gradient towards the lower eastern wall (Figs. 1 and 10).

Conclusion
Cave drip response to unsaturated zone recharge is complex and therefore involves the interaction of several drip pathways with differing response times.This study highlights the importance of hydrogeological controls on water movement in the karst unsaturated zone, which have a critical influence on drip hydrology.The nature of the karst architecture leads to heterogeneous flow in the unsaturated zone, characterized by four different flow types classified using the morphological analysis of Mahmud et al. (2015).This paper applies this method to identify flow types for the individual stalactite discharges measured by continuous hydrological monitoring in a SWWA karst, where hydrological variations are strongly controlled by seasonal variations in recharge.We discover that the discharge data and the morphology-based flow classification agree with each other in terms of flow and geometrical characteristics of ceiling stalactites.We further investigate the drip rates and cave discharge relationship.The mean annual infiltration is found to be 60-70 % of the annual precipitation.Deviations from expected seasonal discharge characteristics have been noted in a few drip discharges.Moreover, the slow dripping sites (icicle and soda-straw) show significant drip variations even while having a uniform nature of drip profiles, indicating differential pressure heads and substantial flow path variability in the overlying unsaturated zone.
LiDAR measurement is used to image the cave ceiling including individual stalactites.Statistical and morphological analyses of the point clouds produced by LiDAR are then used to categorize the ceiling features into different flow types.Later these flow classifications are correlated with the drip logger data to estimate the amount of cave discharge.Moreover, based on the LiDAR coordinates, we paired the locations of individual stalactites with drip loggers on the cave floor.A major advantage of using LiDAR is to upscale/extrapolate the information to the entire cave system to estimate large-scale cave discharge, based on remote sensed data and a minimal number of drip loggers.
We observe no significant relationships between the drip logger mean discharge, skewness, and COV with overburden thickness, due to the possibility of potential unsaturated zone storage volume and increasing complexity of the karst architecture.However, these properties can be used to characterize different flow classes.The correlation matrix shows that similar flow categories are positively correlated, with a significant influence of spatial distribution.We perform the first ever attempt at recharge estimation from cave drip data.Our recharge estimate is quantitative, based on discharge data, and it is in line with measured infiltration from surface rainfall.Therefore, in terms of the mass balance, the method seems to perform well at least at this cave site.
We hypothesize that the amount of discharge from the chamber 1 monitored area is equal to the unsaturated zone recharge within the area of limestone formation.The adjacent chamber ceiling portion is dry around the sites, suggesting that infiltration is being focussed on the studied area of each chamber, which is further agreed on by our conceptual model (Fig. 10b and c).We find that cave discharge per m 2 area is larger than the average surface infiltration (i.e.60-70 % of the rainfall) above the two wettest parts of the chambers and signify water flowing from the surrounding preferential paths to these areas.Infiltration is being focussed to the studied areas of each chamber, by approximately 120 and 170 % respectively.This concentrated flow behaviour is the reason for the contrast with the average surface infiltration values.The majority of the cave ceiling portion is dry at our cave site, suggesting the possibility of capillary effects with water moving around the cave rather than passing through it, especially within studied area of chamber 2 that has relatively lower ceiling elevation and is adjacent to the cave wall.The methodology developed in this paper allows the estimation of deep infiltration without measuring rainfall or ET.Usually, recharge is calculated based on rainfall and ET.Here, we estimate the amount of cave recharge, so the method could be useful to estimate ET, which is difficult to measure.Moreover, this morphology-based flow classification technique can be applied to other cave sites to iden-tify highly focussed areas of recharge, and can help to better estimate the total recharge volume.Nowadays, Zebedee, which is a hand-held LiDAR system (Zlot and Bosse, 2014b, a;Bosse et al., 2012), and drone-mounted systems (Kaul et al., 2016) are available, which make LiDAR technology very useful for three-dimensional mapping.Therefore, Li-DAR can be easily acquired at most of the cave sites if there is possible human access.Moreover, drip loggers have been commercially available for many years (Collister and Mattey, 2008) and are relatively cheap.Therefore, we propose relating these two emerging techniques, so that infiltration can be estimated based on LiDAR data and a minimal number of drip loggers.LiDAR can also be a handy tool in study sites which already have drip data time series (e.g.Jex et al., 2012;Markowska et al., 2015).We demonstrate that morphological properties of stalactites and drip rate monitoring are a suitable means by which to classify karst flow behaviour at a small scale, and should be the focus of future studies using more spatial LiDAR data and temporal drip logger data, limited only by the size of the cave, and a wider range of limestone geologies.This would reduce the uncertainty involved due to the limited scale of this study and would allow for accurate identification of large-scale flow variability.At Golgotha Cave, the collection of drip logger data is ongoing, as well as analysis of tracers of water movement such as stable isotopes.These data will be the focus of future research to expand the possibility of classifying geochemical properties of drip regimes covering large-scale observation, and even describing the hydrodynamic response of the unsaturated zone in the cave (Carrasco et al., 2006;Blondel et al., 2010;Mudry et al., 2008).

Figure 2 .
Figure 2. (a) Box plot of monthly rainfall at the Golgotha Cave site.(b) Monthly water budget for 3 observed hydrological years.

Figure 3 .
Figure 3.Both monitored ceiling images from T-LiDAR with 34 identified stalactites feeding drip loggers, (a) chamber 1 and (b) chamber 2. Blue circles indicate locations of the ceiling under which the loggers are placed (indicated by red arrows).The blue arrows in both figures show the geographic orientation.These ceiling stalactite locations are identified based on both on-site field observation and the alignment analysis of drip loggers and stalactites (Sect.4.2).

Figure 4 .
Figure 4. Morphological analysis of the cave ceiling portion consisting of sites 1A, 1B, and 1i used to identify different flow patterns (icicle, soda straw, fracture, and combined).(a) Cave ceiling topography in 2-D.Colour scales represent elevations in metres relative to the T-LiDAR receiver.Stalactites identified in (b) are all categorized as icicle flow pattern(Mahmud et al., 2015).

Figure 5 .
Figure 5. Morphological analysis of the cave ceiling portion consisting of sites 1ii, 1iii, and 1v to identify different flow patterns (icicle, soda straw, fracture, and combined).(a) Cave ceiling topography in 2-D.Colour scales represent elevations in metres relative to the T-LiDAR receiver.(b) Stalactites categorized as icicle flow pattern.(b) Stalactite 1v categorized as soda straw(Mahmud et al., 2015).

Figure 6 .
Figure 6.Both chambers' drip rate time series for the entire monitoring period.(a) Chamber 1 drip rates.Further classification of chamber 2 drip sites for effective visualization: (b) slow flow rates with drip frequency of less than 10 per 15 min, (c) medium discharges with drip frequency between 10 and 100 per 15 min, and (d) fast drip rates with more than 100 drips per 15 min.

Figure 7 .
Figure 7. Actively dripping stalactite count.Green circles are stalactites identified as active in terms of dripping (displaying water droplets and/or white flash at their tip) and yellow ones are inactive.In this photo, 32 green circles and 13 yellow ones are counted.

Figure 8 .
Figure 8. Drip data characteristics for different flow classifications.

Figure 9 .
Figure 9. Correlation matrix for various drip sites with flow type classification for (a) chamber 1 and (b) chamber 2. Different flow types (I: icicle; F: fracture; C: combined; S: soda straw) characterized by the morphological analysis are shown in parentheses (left y axes and x axes).Limestone thicknesses over the drip sites are on the right y axes.The colour scale indicates the correlation coefficient between drip time series.

Figure 10 .
Figure 10.Representation of groundwater flow for the Tamala limestone formation in Golgotha Cave.(a) shows a vertical section of the entire cave with possible groundwater flows.(b) and (c) illustrate close views of chambers 1 and 2 labelling areas of ground surface infiltration and drip discharge from cave ceilings (in red circle).The "?" and arrows pointing to the left and right in (b) and (c) show that we cannot be sure of the precise flow route, although we are able to quantify the flow amount.

Table 1 .
Flow classification of chamber 1 drip data.

Table 2 .
Flow classification of chamber 2 drip data.

Table 3 .
Cave discharge calculation for both chambers.