Groundwater seepage is a key driver of theater-headed valley formation in limestone

Whereas field observations and experimental/numerical simulations have shown that groundwater seepage can lead to the formation of theater-headed valleys (THVs) in sand and gravel, the mechanism for their development in consolidated bedrock remains controversial (Lamb et al., 2006; Lapotre and Lamb, 2018; Lobkovsky et al., 2007; Micallef et al., 2021; Ryan and Whipple, 2020). The classic model for the formation of THVs in bedrock entails weathering at the seepage face, which weakens the rock and renders grains cohesionless; these grains are subsequently removed by seeping water, leading to undermining, collapse, and retreat (Dunne, 1980; Laity and Malin, 1985). Arguments against the role of groundwater in the formation of THVs in clastic sedimentary or non-sedimentary rocks include: (1) the large size of some THVs, (2) time required for extensive weathering, and (3) inability of spring discharge to remove the talus (Lamb et al., 2006; Lapotre and Lamb, 2018). Groundwater weathering and erosion processes in different types of bedrock have not been quantified in terms of mechanisms and rates. A diagnostic link between groundwater seepage and geometry of bedrock THVs remains unclear and poorly quantified (Petroff et al., 2011). Reasons for such gaps in knowledge include: (1) field observations and experimental simulations are difficult—the time scales involved are long and the processes are complex; (2) experimental and numerical analyses have been based on simplistic assumptions about flow processes and hydraulic characteristics, and they rarely simulate effects of geologic heterogeneities; and (3) an unambiguous example of a THV formed by groundwater seepage in bedrock is lacking because of overlapping features generated by rainfall runoff. In addition, bedrock THVs are not uniquely associated with groundwater. Groundwater seepage may be complemented by surface-water erosion, or the latter may be enough to explain bedrock THV formation, e.g., in substrate with stratigraphic contrasts and vertical joints, waterfall erosion through plunge-pool erosion or plucking/toppling, and megaflooding (Lamb et al., 2008, 2007, 2006, 2014; Nash, 1996; Ryan and Whipple, 2020; Scheingross and Lamb, 2017). Elucidating the association between groundwater seepage and bedrock THVs is crucial to reconstructing/predicting landscape evolution, with direct relevance to hazard assessment, and to inferring past climates and fluvial conditions.

In this study, we combine field and remote-sensing observations from the Maltese Islands with numerical modeling to demonstrate that groundwater seepage can be the main driver of THV formation in limestone. We infer the mechanisms responsible and the factors controlling THV location, shape, and size.

The Maltese Islands comprise four main Oligo-Miocene sedimentary formations (Pedley et al., 1976; Fig. 1A). Lower Coralline Limestone (Chattian) and Globigerina Limestone (Aquitanian to Langhian) consist of shallow-water marine calcareous algae and deep-water planktonic foraminifera, respectively. Blue Clay (Langhian to Tortonian) is composed of banded kaolinitic marls and clays. Upper Coralline Limestone (late Tortonian to early Messinian) is subdivided into four members: Gebel Imbark, Tal-Pitkal, Mtarfa, and Ghajn Melel. Tal-Pitkal Member is a coarse-grained wackestone and packstone containing coralline algal, mollusc, and echinoid bioclasts. Mtarfa Member consists of massive to thickly bedded carbonate mudstone and wackestone rich in coralline algae. The entire sedimentary sequence is disrupted by two normal fault systems, the most widespread of which consists of ENE-WSW–trending faults dipping at 55°–75° (Illies, 1981).

The Maltese Islands have a semiarid climate, with a mean annual precipitation of 550 mm (Galdies, 2011). Lower Coralline Limestone hosts the mean-sea-level groundwater, whereas Upper Coralline Limestone hosts perched groundwater bodies above Blue Clay (Bakalowicz and Mangion, 2003).

All THVs considered in this study occur in Upper Coralline Limestone (Fig. 1A; Haroon et al., 2021). Onshore THVs are associated with natural springs fed by perched groundwater, and most of them are oriented parallel to known faults. One of these THVs—Gnejna Valley in northwestern Malta (Fig. 1A)—was selected as a case study because of its size and accessibility. The mean temperature, conductivity, and pH of perched groundwater at Gnejna Valley for 2009–2017 CE are 19.7 °C, 1420 μS/cm, and 7.5, respectively. Offshore THVs are thought to have formed subaerially because they were exposed during the Last Glacial Maximum and are draped by 5–10 m of hemipelagic sediment accumulated since then (Micallef et al., 2013).

Topographic lidar and multibeam echosounder data are available for the entire Maltese Islands and for the seafloor to as much as 150 m depth, respectively (Micallef et al., 2013, 2019; Fig. 1A). These data were used to map onshore and offshore THVs. Gnejna Valley was surveyed with an uncrewed aerial vehicle to generate orthophotos and elevation models with a horizontal resolution of 10 cm/pixel (Fig. 1B; Fig. S1 in the Supplemental Material¹). Electrical resistivity tomography (ERT) and ground penetrating radar (GPR) data were acquired along two 470-m-long profiles located to the north and east of Gnejna Valley, respectively, to map the underlying stratigraphy, structure, and groundwater (Fig. 1B). A 250-m-long scanline survey was carried out at the base of the head of Gnejna Valley (Fig. 1B) to estimate orientation, spacing, and roughness of discontinuities (ISRM, 1978, 1985). Rock-wall strength was measured with an N-type Schmidt hammer. Surface samples of Upper Coralline Limestone and Blue Clay were taken from the head and northern rim of the valley for scanning electron microscopy, measurement of geotechnical/hydraulic properties, and cosmogenic nuclide dating (Fig. 1B).

Hydrological modeling entailed estimation of (1) groundwater discharge at Gnejna Valley using GIS hydrology tools and (2) surface discharge required to transport boulders in the valley using continuity, Manning, and Shields equations. A landscape evolution model, based on Chen et al. (2014), was developed to assess the effect of erosion by surface flow on valley development. Three-dimensional (3-D) distinct element modeling was carried out on a 1000 m × 500 m × 40 m domain using 3DEC modeling code (Fig. 1C) and accounted for elastic hydromechanical widening, dissolution of joints/fractures in limestone, sliding of limestone blocks, and creep in clay (Barton, 1972; Garcia-Rios et al., 2015; Hoek and Bray, 1981; Roscoe and Burland, 1968).

Detailed information on the data and methods used is available in the Supplemental Material.

We identified 53 THVs: 30 across Malta and Gozo, and 23 offshore the northeastern coast (Fig. 1A). THVs have U-shaped cross-sections and linear longitudinal profiles; valleys widen downstream at a mean rate of 0.67 m width per 1 m length (σ = 0.65). With a length of 1.5 km, width of 0.85 km, and relief of 75 m, Gnejna Valley is the largest THV (Fig. 1B). The valley head and flanks are formed in Tal-Pitkal and Mtarfa Members dipping at <1° to the west-northwest, whereas the bed predominantly consists of Blue Clay (Fig. 2A). The head of Gnejna Valley is a near-vertical, 25-m-high wall characterized by vertical joints and fractures (Fig. 2A; Fig. S3; geomechanical properties in Table S2 in the Supplemental Material). At the base of the head, water seeps out of joints and fractures and a talus of Upper Coralline Limestone boulders, mainly >1.5 m in diameter, has formed (Fig. 2A). Mtarfa Member outcropping at the valley head and flanks is characterized by a recess with evidence of spalling (Figs. 2B and 2D). Dissolution occurs only in the inner faces of joints and fractures (Fig. 2C). Surface exposure ages for samples 01 and 02 are 11.1 and 33.4 ka, respectively (Table S6).

Electrical resistivity tomography (ERT) profiles show a shallow, 10–30 m high-resistivity layer (>600 Ωm), a ~5 m intermediate-resistivity layer (100–300 Ωm), and a deeper, low-resistivity layer (<50 Ωm) (Fig. 2G; Fig. S2). Via correlation with outcrops, we interpret these layers as unsaturated Upper Coralline Limestone, saturated Upper Coralline Limestone (i.e., containing perched groundwater), and Blue Clay, respectively. The southeastern section of profile T4 shows an offset that corresponds with a west-east–trending fault throwing 15 m toward the south (Fig. 2G) and with the location of the main spring. GPR profiles show quasi-vertical master joints and fractures in Upper Coralline Limestone and a decrease in electromagnetic signal amplitude associated with perched groundwater (Fig. 2G; Fig. S2).

Based on the above data, we dismiss three hypotheses for THV formation at Gnejna Valley:

1. Detachment of valley head by surface flow: Apart from a 20-m-wide and 90-m-long V-shaped channel above the valley head (Fig. 2E), there is a lack of well-developed channels draining toward Gnejna Valley (Fig. 1A; Fig. S1). Fluvial abrasion marks and plunge pool undercuttings are absent (Figs. 2E and 2F). Landscape evolution modeling suggests that erosion by surface flow would form a valley with a pointed head, which is different from Gnejna Valley (Fig. S5).

2. Erosion by megafloods: Features cited in support of megaflooding, such as scours and flood deposits, are absent. In addition, the orientation of THV axes is variable (Fig. 1A). There are no records of megafloods across the Maltese Islands, and factors responsible for such events elsewhere (e.g., rapid melting of ice sheets, ice-dammed lake failure) have been absent here in the last 5 m.y.

3. Groundwater weathering at a seepage face, leading to undermining and loss of support: Groundwater seeps through the joints/fractures rather than the low-permeability matrix. Recesses are pervasive across Mtarfa Member outcrops everywhere, including other sites and not just at the head of Gnejna Valley, and their formation is likely associated with mechanical weathering by wetting and drying (Figs. 2B and 2D).

We propose an alternative mechanism by which Gnejna Valley was formed (Fig. 3). Development and upslope retreat of the THV head in Upper Coralline Limestone takes place via slope failure—toppling or sliding—driven by groundwater seepage via joints and fractures. The latter gives rise to: (1) widening of valley head–orthogonal joints and fractures by fluid pressure and dissolution, which reduces Upper Coralline Limestone's mass strength; and/or (2) creep of Blue Clay, which leads to loss of support at the base of the THV head. Failed material is transported downslope and away from the head by creep of Blue Clay (driven by seepage and precipitation) and sliding. We test the effectiveness of this mechanism in forming THVs using 3-D distinct element modeling.

We consider the above two groundwater-related processes individually and together. Given that faulting appears to control the location of seepage (Fig. 2), we also assess the role played by fault-zone permeability (Caine et al., 1996). Four modes of groundwater seepage are examined: (1) point seepage (at the fault), and (2) no change, (3) linear decrease, and (4) exponential decrease in seepage away from the fault and across the fault damage zone (Fig. 1C). The latter mode is based on the exponential increase in fracture spacing with fault-normal distance, which is commonly reported in the literature and observed at Gnejna Valley (Liao et al., 2019; Sagy et al., 2001; Fig. S3). Simulations taking into account increased seepage across the fault damage zone (modes 2–4) result in valleys similar to THVs across the Maltese Islands, widening downstream at rates of 0.6–2 m width per 1 m length (Fig. 4). However, simulations for mode 4 give rise to a head morphology that is most similar to that of Gnejna Valley and most other THVs. The width of the fault damage zone controls the maximum width of the THV. Both widening of joints and fractures in Upper Coralline Limestone and creep of Blue Clay are effective processes driving THV formation. The estimated THV head retreat rate is 0.72–1.32 cm/yr. Surface exposure ages suggest retreat rates of 1.8–7.1 cm/yr (calculated by dividing axial distance of sampling sites from THV head by age). We would expect actual retreat rates to be closer to 1.8 cm/yr because the young age recorded at sampling site 01 was likely due to a recent slope failure.

The talus accumulating on the THV bed needs to be removed to ensure that head retreat is not hindered (e.g., Dunne, 1990; Lapotre and Lamb, 2018). The 3-D distinct element model yields creep and sliding rates of 1–62 cm/yr for a 1.5 m Upper Coralline Limestone boulder on Blue Clay sloping 1°–20° (Fig. S6). Considering that the slope gradient of Blue Clay in Gnejna Valley is 2°–6.3° and that it likely reached 20° at the coast prior to THV formation, the transport rates of 1.5 m boulders on the Gnejna Valley bed have generally been higher than head retreat rates. Spring discharge (estimated at 0.01 m³/s) is too low to remove these boulders. A flash flood generated by a rainfall event of at least 10.3 m/day would be required to transport such boulders (Supplemental Material section 3.1), which is unrealistic. The largest boulders that may have been transported by documented flash floods (rainfall of 100 mm/day) are ~10 cm in diameter.

Estimates of rates for THV head retreat and boulder transport are based on current annual precipitation. Past climate variability in the Maltese Islands is poorly constrained, and its exclusion limits the accuracy of our estimates. On the other hand, factors promoting head retreat (e.g., increased annual rainfall) are expected to promote removal of talus as well.

We demonstrate how groundwater seepage can be the main driver of THV formation in jointed limestone overlying clays. By integrating field and remote-sensing observations from the Maltese Islands with numerical modeling, we show that groundwater seepage due to a permeability contrast can widen joints and fractures in limestone and promote creep of the underlying clay. Both of these processes lead to development and retreat of a THV head by slope failure. Talus is removed by creep and sliding on the valley bed at a rate faster than that of head retreat. By concentrating groundwater seepage across a cliff face, faults control the location of THVs. The maximum width of the THV is determined by the width of the fault damage zone, whereas the shape of the head is controlled by how seepage changes away from the fault core. An exponential decrease of seepage gives rise to a theater-shaped head that is similar to that of Gnejna Valley. Our model may explain the formation of THVs by groundwater in jointed, strong-over-weak chemical sedimentary lithologies, particularly in arid terrestrial settings.

This project was funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement 677898 [MARCAN]) and HYFREW (CNR–University of Malta). We thank the three anonymous reviewers for their constructive comments.