Box canyon erosion along the Canterbury coast (New Zealand): A rapid and episodic process controlled by rainfall intensity and substrate variability

Box canyon formation has been associated to groundwater seepage in unconsolidated sand to gravel sized sediments. Our understanding of box canyon evolution mostly relies on experiments and numerical simulations, and these rarely take into consideration contrasts in lithology and permeability. In addition, processbased observations and detailed instrumental analyses are rare. As a result, we have a poor understanding of the 20 temporal scale of box canyon formation and the influence of geological heterogeneity on their formation. We address these issues along the Canterbury coast of the South Island (New Zealand) by integrating field observations, optically stimulated luminescence dating, multi-temporal Unmanned Aerial Vehicle and satellite data, time-domain electromagnetic data, and slope stability and landscape evolution modelling. We show that box canyon formation is a key process shaping the sandy gravel cliffs of the Canterbury coastline. It is an episodic 25 process associated to groundwater flow that occurs once every 227 days on average, when rainfall intensities exceed 40 mm per day. The majority of the box canyons in a study area SE of Ashburton has undergone erosion, predominantly by elongation, during the last 11 years, with the most recent episode occurring 3 years ago. The two largest box canyons have not been eroded in the last 2 ka, however. Canyons can form at rates of up to 30 m per day via two processes: the formation of alcoves and tunnels by groundwater seepage, followed by retrogressive 30 slope failure due to undermining and a decrease in shear strength driven by excess pore pressure development. The location of box canyons is determined by the occurrence of hydraulically-conductive zones, such as relict braided river channels and possibly tunnels, and of sand lenses exposed across sandy gravel cliff. We also show that box canyon formation is best represented by a linear diffusive model and geometrical scaling. 35 https://doi.org/10.5194/esurf-2020-29 Preprint. Discussion started: 2 June 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Whereas erosion by overland flow has been extensively studied and associated geomorphic rate laws have been developed (Sklar and Dietrich, 2001;Whipple et al., 2000), erosion by groundwater flow has received considerably less attention. Groundwater can mechanically erode or deform sediment and bedrock by removal of mass from a seepage face through flowage, frictional slipping and particle dislodgement (Dunne, 1990;Iverson and Major, 1986). Groundwater seepage can also lead to slope instability by undercutting, reduction of soil shear strength and shear surface development (Carey and Petley, 2014;Chu-Agor et al., 2008). Groundwater has been 45 implicated as an important geomorphic agent in valley network development, both on Earth and on Mars (Abotalib et al., 2016;Dunne, 1990;Harrison and Grimm, 2005;Higgins, 1984;Kochel and Piper, 1986;Malin and Carr, 1999;Salese et al., 2019). For over one hundred years, a characteristic morphology -comprising theatre-shaped heads, steep and high valley walls, constant valley width, flat floors, low drainage densities, and short tributaries with large junction angles -has been cited as diagnostic of groundwater activity (Abrams et al., 2009;Dunne, 50 1990;Higgins, 1984;Russel, 1902;Schumm and Phillips, 1986). We refer to valleys with these characteristics as box canyons.
The classic model explaining the development of box canyons entails a canyon headwall that lowers the local hydraulic head and focuses groundwater flow to a seepage face. This leads to upstream erosion by undercutting, 55 the rate of which is limited by the capacity of seepage water to transport sediment from the seepage face (Abrams et al., 2009;Dunne, 1990;Howard and McLane, 1988). Groundwater seepage has been shown to unambiguously lead to box canyon formation in unconsolidated sand to gravel sized sediments (Dunne, 1990;Lapotre and Lamb, 2018) (Table 1). In sediments finer than sands, erosion is typically limited by detachment of the grains at the seepage face. In silts and clays, the permeability is so low that the groundwater discharge is often less than that required to overcome the cohesive forces of the grains (Dunne, 1990). The role of groundwater seepage and box canyon formation in bedrock, on the other hand, remains controversial (Lamb et al., 2006;Pelletier and Baker, 2011). Gravel braided river deposits (Schumm and Phillips, 1986) Florida Panhandle, USA Non-marine quartz sands that contain discontinuous layers of clay or gravel (Schumm et al., 1995) Kalahari, southern Africa Conglomerates, marls, duricrusts and unconsolidated sands (Nash et al., 1994) Martha's Vineyard and Nantucket Island, USA Outwash gravelly sand (Uchupi and Oldale, 1994) Obara area, Japan Granodiorite regolith (Onda, 1994) South Taranaki, New Zealand Terrestrial dune sand and tephra overlying marine sands and gravels (Pillans, 1985) Vocorocas, Brazil Alluvial sands (Coelho Netto et al., 1988) Our understanding of box canyon evolution is predominantly derived from theoretical, experimental and numerical models, which have reproduced the characteristic box canyon morphologies in unconsolidated sediments (Chu-Agor et al., 2008;Higgins, 1982;Howard, 1995;Lobkovsky et al., 2004;Pelletier and Baker, 70 2011;Petroff et al., 2011;Wilson et al., 2007). Such studies suggest that the velocity at which box canyon heads advance is a function of groundwater flux and the capacity of seepage water to transport sediment from the seepage face (Abrams et al., 2009;Fox et al., 2006;Howard and McLane, 1988), and that canyon head erosion occurs by episodic headwall slumping (Howard, 1990;Kochel et al., 1985). Streams incised by groundwater seepage have been shown to branch at a characteristic angle of 72° at stream tips, which increases to 75 120° near stream junctions Yi et al., 2017), whereas growing indentations competing for draining groundwater result in periodically-spaced valleys (Dunne, 1990;Schorghofer et al., 2004). Canyon network geometry appears to be determined by external groundwater flow field rather than flow within the canyons themselves .

80
A number of fundamental questions related to the evolution of box canyons in unconsolidated sediments remain unanswered. Firstly, the temporal scale at which box canyons form is poorly quantified due to a paucity of processbased observations and detailed instrumental analysis. Field observations of groundwater processes are rare (e.g. Onda, 1994), primarily due to the difficulty with accessing the headwalls of active canyons, the potential long timescales involved, and the complexity of the erosive process (Chu-Agor et al., 2008;Dunne, 1990). Quantitative 85 assessments of box canyon evolution have relied on experimental and numerical analyses, but these tend to be based on simplistic assumptions about flow processes and hydraulic characteristics. Experimental approaches are based on a range of different methods, which limits comparison of their outcomes (Nash, 1996). Published erosion rates vary between 2-5 cm per century (Abrams et al., 2009;Schumm et al., 1995) and 450-1600 m 3 per year (Coelho Netto et al., 1988). Secondly, the influence of geologic heterogeneities on box canyon evolution is also 90 poorly understood. Lithological strength and permeability contrasts are rarely simulated by experimental and numerical analyses. Thirdly, there only a few places where the mechanisms by which seepage erosion occurs have been clearly defined (e.g. Abrams et al., 2009). Basic observations and measurements of box canyon erosion rates and substrate geologic heterogeneities are needed to test and quantify models for canyon formation and improve our ability to reconstruct and predict landscape evolution by groundwater-related processes.

95
In this study we revisited the Canterbury Plains study site (Schumm and Phillips, 1986) and carried out field observations, geochronological analyses, repeated remote sensing surveys, near-surface geophysical surveying and numerical modelling of coastal box canyons to: (i) identify the processes by which groundwater erodes box canyons along the coast, (ii) quantify the timing of box canyon erosion and its key controls, and (iii) assess the 100 influence of geological/permeability heterogeneity on the box canyon formation process.

Regional setting
The flat to gently inclined Canterbury Plains, located in the eastern South Island of New Zealand, extend from sea 105 level up to 400 m above sea level, and cover an area 185 km long and 75 km wide (Fig. 1a). A series of high energy braided rivers emerge from the >3500 m high Southern Alps and flow south-eastwards to the shoreline (Kirk, 1991). The plains were formed by coalescence of several alluvial fans sourced from the these rivers (Browne and Naish, 2003;Leckie, 2003). The Quaternary sedimentary sequence comprises a >600 m thick succession of cyclically stacked fluvio-deltaic gravel, sand and mud with associated aeolian deposits and 110 palaeosols (Bal, 1996;Browne and Naish, 2003). The gravels consist of greywacke and represent a variety of channel fill beds and bar forms, whereas the isolated bodies of sand are relict bars and abandoned channels. The interglacial sediments are better sorted and have higher permeability than the glacial outwash, resulting in a wide range of hydraulic conductivities (Scott, 1980). New Zealand's largest groundwater resource is hosted in the gravels down to at least 150 m depth (Davey, 2006). The upper Quaternary sediments are exposed along a 70 km 115 long coastline south-west of the Banks Peninsula (Moreton et al., 2002). This coastline is retrograding at approximately 0.5-1 m per year and consists of cliffs fringed by mixed gravel and sand beaches (Gibb, 1978).
The study area is a 2.5 km long stretch of cultivated coastline located 16 km to the south-east of Ashburton (Fig.   1b). The coastline within the study area consists of a 15-20 m thick exposure of poorly-sorted and uncemented matrix-supported outwash gravel, which is capped by up to 1 m of post-glacial loess and modern soil (Berger et 120 al., 1996). The cliff face is punctuated by ~0.5 m thick lenses of sand or clean gravel.   Fig. 1b). The latter were collected by hammering stainless steel tubes into the sediment and ensuring that the material was not exposed to light. 150 location determined by differential GPS. Orthophotos and digital elevation models with a horizontal resolution of 10 cm/pixel were generated from the UAV data using Drone Deploy. The surveys were carried out on the   and Macnae (1991) and Fitterman (2015). The survey parameters included 4 turns, a 10 × 10 m 2 square TX loop, 160 and a TX current output of 1 A. The G-TEM was operated in a fixed offset-sounding configuration, which is termed "Slingram" mode, in which the RX coil was placed 30 m from the centre of the TX loop and the TX-RX pair moved together along a linear transect at 5 m station spacing, maintaining the 30 m offset. All soundings were collected in the 20-gate mode with an acquisition interval of 6 × 10 -6 s to 8 × 10 -4 s (after ramp-off), corresponding to investigation depths of ~80 m. At each station, a consistent 1-D smooth model of electrical 165 resistivity vs. depth was performed based on the iterative Occam-regularised inversion method (Constable et al., 1987) and using the IXG-TEM software from Interpex Limited.

170
Satellite images with a horizontal resolution of 1 m/pixel and dating back to 2004 were obtained from Google Earth. Precipitation records dating back to 1927 were provided by Environment Canterbury. The latter also provided a time series of water level data since 2015 from a 30 m deep well located 10 km to the north-east of the study area.

Sample analyses
Sediment samples were analysed for grain size distribution using sieves following the ASTM D0422 protocol.

180
The composition of the coating on selected sediment outcrops within the box canyons was determined using Xray fluorescence.

190
In order to obtain luminescence ages, two types of measurements were performed. The dose accrued by the crystal from natural radioactivity since its last exposure to daylight (called the palaeodose) was determined as an equivalent dose (De). This was done by measuring the light emission of the crystal upon optical stimulation, and matching this emission to signals generated by the exposure to a known dose of radiation given in the laboratory.

195
This is expressed as the amount of absorbed energy per mass of mineral (1 J kg -1 = 1 Gy (Gray)). Radioactivity measurements were carried out on each sample in order to determine the annual dose (Da), which represents the rate at which the environmental dose was delivered to the sample (Gy ka -1 ). The age was obtained by dividing the two determined parameters. As low luminescence sensitivity and poor dosimetric characteristics were reported for quartz from sediments in New Zealand (see Preusser et al. (2009)

Morphological change detection
The method used to measure box canyon aerial erosion in between surveys entailed the manual delineation of shapefiles around box canyon boundaries for each survey (using orthophotos, digital elevation models and slope gradient maps in the case of the UAV data, and satellite images in the case of the Google Earth data), the estimation of their areas, and the comparison of the latter in between surveys. The uncertainty inherent in this approach is 215 related to the digitisation of the canyon boundaries. We made sure that a vertex was added at least every 5 pixels for both the UAV and Google Earth data. This ensures that a minimum erosion of 0.25 m 2 (in the case of the UAV data) and 25 m 2 (in the case of the Google Earth data) was detected.

220
Slope stability modelling We developed a slope stability model based on the limit equilibrium and segmentation strategy of the Bishop method, where a soil mass is discretised into vertical slices. The factor of safety Ff is calculated using the following 225 (Fredlund and Krahn, 1977;Fredlund et al., 1981): where c' (in kPa) is effective cohesion, ɸ' (in °) is effective angle of friction, u (in kPa) is pore-water pressure, D

230
(in kN) is concentrated point load, β (in m) represents geometric parameters, ω is geometric parameter, and α (in °) is inclination of slice base. N is the normal force acting on the slide base and can be computed by: where W (in kN) is slice weight (unit weight γs (in kN m -3 ) × volume (in m 3 )) and k is hydraulic conductivity (in m s -1 ).
We also modelled the water flow and pore pressure distribution within the soil using the Poisson equation, which is the generalised form of the famous Laplace equation (Whitaker, 1986): where q is the total discharge (m 3 s -1 ), kx and ky are equal to the hydraulic conductivity (m s -1 ) in the horizontal and vertical directions, respectively, and h is the hydraulic head (m).
Equation (3) applies to water flow under steady-state and homogeneous conditions, whereas the following equation is applicable to dynamic and inhomogeneous conditions: where describes how the volumetric water content changes over the time.
The water transfer theory accounts for transient behaviour, which can be defined by the following equation (Domenico and Schwartz, 1997): where ̇ is the cumulative mass of water that enters the porous medium, ̇ is equal to the mass of water that leaves the porous medium, and Ms is the mass source within the representative elementary volume. The rate of 260 increase in the mass of water stored within the representative elementary volume is: where ̇ and ̇ represent the rate of change of liquid water and water vapour, respectively.

265
The mechanical and hydraulic soil properties employed in this model are listed in Table 2 and were obtained from Dann et al. (2009) andAqualinc Research Limited (2007). We modelled two scenarios, based on the available rainfall data (see Sect. 4.4). The first is a 3-day long intense rainfall event ((I-D)3) covering the period 20 th July -22 nd July 2017. The second is a 14-day period with occasional, low intensity rain ((I-D)14) between the 21 st June 270 and 4 th July 2017. Each scenario is modelled for two sandy gravel slopes, one with a 0.5 m thick sand lens and the other with a 0.5 m thick gravel lens. Lateral water inflow and surface water infiltration were estimated from the hydrological model in Micallef et al. (2020). Slope stability modelling was carried out using the Slide2 software package by Rocscience.

Landscape evolution modelling
We also built a landscape evolution model (LEM) using the Python modelling environment Landlab (Barnhardt 280 et al., 2020;Hobley et al., 2017). The model allowed us to simulate two main processes: groundwater flow and associated erosion (Fig. 2).

295
We used the model to test two erosional modes. The first mode is the stream power law model, which simulates surface erosion (Barnhart et al., 2019;Braun and Willett, 2013): S is the slope (dimensionless) and m and n are exponents (dimensionless). We assumed m = 0.5 and n = 1, due to a lack of data and a robust methodology to determine these two parameters (Harel et al., 2016).

305
The second mode is a linear diffusion to simulate erosion by gravitational sediment movement (Barnhart et al., 2019;Culling, 1963): For the sake of simplicity and to easily compare the results, K (Eq. (7)) is computed as D (Eq. (8)), which we consider to be proportional to surface water shear stress and seepage flux: where τb is the shear stress [MLT -2 ], τt is the threshold shear stress [MLT -2 ] (assumed to have a value of 0.03 N m -2 ), and M is an empirical parameter [L 3 T -1 ] that we assume to have a value of 0.1 m 3 s -1 . The values for τt and 320 M are empirical and were estimated by trial and error.
Additionally, as a theoretical exercise, we tested the evolution of the surface using the stream power law approach with surface flow and without groundwater seepage. Equation (7) was changed by including surface water discharge (Q) from an upstream drainage area: 12 The initial topography was extracted from a digital elevation model of part of the study area that has not been affected by canyon erosion (Fig. 2b), which allowed us to incorporate a realistic approximation of the cliff morphology. To simulate the groundwater flow, we arbitrarily defined the base of the aquifer at 0.1 m above sea level, to minimise vertical groundwater movement, and applied the Dupuit-Forchheimer approximation. The 335 initial water table was placed at 5.9 m above sea level (in accordance with well data). Since Landlab is a 2-D modelling environment, it was not possible to include sand lenses. For this reason, we simplified the geology by considering two types of sediment: sandy gravels and with a NW-SE strip of sand. Their distribution and properties are shown in Fig. 2c and Table 2, respectively. In Eq. (10), we assumed that the erosion coefficient (K) of the sand unit is the double that of sandy gravel.

340
The aquifer recharge was modelled using two different precipitation scenarios. The 'high intensity' scenario is based on a first hour of low recharge (1 × 10 -9 m s -1 ) and a 2 hour long high recharge storm (2.8 × 10 -5 m s -1 ). The latter value was selected to replicate the rise in the water table reported after the 21 st July 2017 storm (see Sect. 4.4.1). The 'low intensity' scenario consists of a first hour of low recharge (1 × 10 -9 m s -1 ), followed by a 5 hour 345 long low recharge storm (1.12 × 10 -5 m s -1 ). The total modelling time was 24 hours.

Box canyons along the Canterbury coast -distribution and morphology
We have mapped 315 box canyons (locally also known as "dongas") along 70 km of the Canterbury coastline (mean of 4.5 canyons per km of coastline). The distribution of the box canyons is clustered (nearest neighbour ratio of 0.33 with a z-score of -22.67 and a p-value of 0); the majority of the canyons are located between Rakaia and Rangitata Rivers (Fig. 1a), particularly in the vicinity of Ashburton River. The heads of many box canyons 355 connect to shallow, relict meandering channels (Fig. 3a). Some of these channels are visible in aerial photographs, in spite of the terrain having been worked by farmers (Fig. 3b). a length to width ratio that varies between 1 and 7.9, with a mean of 2 (standard deviation of 0.89) (Fig. 3c).

Field site observations
In May 2017, our study area hosted 33 box canyons that vary between 15 m and 600 m in length (Figs. 1b; 4a).

375
During the site visits we did not encounter evidence of surface flow. However, the middle to lower sections of the canyon walls and cliffs were consistently wet. These sections were also characterised by failure scars and alcoves, particularly above the sandy lenses. Alcoves were also encountered at the base of canyon heads, where they were wet and up to 1 m deep (Fig. 4e). Some sandy layers outcropping across the cliff face hosted tunnels ( Fig. 4c-d).
Above these tunnels, theatre-shaped scars with a shallow and narrow gully at their base were observed (Fig. 4c).

380
At the base of the scars, the box canyon heads and some box canyon mouths, we encountered mass movement debris that was occasionally intact and that predominantly consisted of gravel, sandy gravel and loess (Fig. 4c, h).
Box canyons have gravel-covered irregular floors. Whereas the smaller box canyon have a U-shaped crosssection, the larger box canyons have gently sloping V-shaped cross-sections, with loess draping their walls (Fig.   4f). Sandy and clean gravel layers outcropping within the canyons were wet; the former appeared weathered, whereas the latter were coated by Fe and Mn (Fig. 4g). Fences were locally seen suspended across a number of box canyons (Fig. 4b).

Luminescence ages
Four sets of OSL ages are presented in Table 3. The pIRIR290 ages are higher than the ages obtained by applying pIRIR225 protocol. The cause of this difference is not yet fully understood, although it can partially be attributed to the results of the dose recovery test and the poor bleachability of the pIRIR290 signals compared to pIRIR225 signals . Considering that no anomalous behaviour of the investigated signals was observed (see Supplementary Materials), we are unable to explain the overestimation of the K-feldspar ages compared to the polymineral fine grain ages in the case of NZ13A, especially since the opposite behaviour is observed in the 400 case of sample NZ14A. However, considering a 95% confidence level, ages obtained using different methods broadly overlap, the only exception being the pIRIR225 ages obtained on K-feldspars on sample NZ14A, which we regard as an outlier.   Figure 6 shows the total area eroded between surveys (which amounts to 3273 m 2 ), the daily precipitation and the associated changes in water table height. Only three surveys recorded box canyon erosion. Two of these surveys happened soon after rainfall events of >40 mm in one day (Fig. 6). The most important of these covers the period between the 15 th and 23 rd July

480
The first-time-gate profile of transect May15-1 is located upslope of small but recently eroded box canyons (Fig.   9a). Near the middle of this transect there is a distinctive peak that is much higher than the background. The peak is ~20-30 m wide and it appears in a similar fashion on each of the gates 1 through 7 (not shown here), although it cannot be clearly observed after gate 7. Transect May 15-2 is located upslope of recently eroded box canyons in the south-west and relatively less active canyons in the north-east of the investigated area (Fig. 9b). Lateral comparison to the previous two transects (Fig. 9c). Based on all three profiles, a general observation that can be made is that the first-time-gate amplitude of the slingram response is higher upslope of the more recently active box canyons. The factory of safety of the slope prior to any rainfall event was 2.514. During the first scenario ((I-D)3), the factor of safety decreased to 1.371 due to undermining by tunnelling associated to high pore pressures within the sand lens, and then to 0.614 as a result of a decrease in the shear strength of the lower slope material due to an increase 510 in pore pressure (Figs. 10a, 11a). A rainfall intensity of 40 mm per day is required to bring the factor of safety below 1 (Fig. 11c). In the case of the second scenario ((I-D)14), changes in pore-water pressure did not result in either tunnelling or slope failure. This only resulted in a decrease in the effective stress and in the factor of safety (1.216) (Figs. 10b, 11b).  The factor of safety of the slope prior to any rainfall event is 1.793. For the first scenario ((I-D)3), the factor of safety decreased to 1.166, and neither tunnelling nor slope failure occurred (Figs. 12a-b). In the case of the second scenario ((I-D)14), the outcome is the same, with the factor of safety decreasing to just 1.588 (Figs. 12c-d).

LEM
The model results for the two erosional modes and two rainfall scenarios are displayed in Figs. 13-16 and 540 discussed in more detail below.
The results of the stream power law model with a high intensity rainfall scenario are shown in Fig. 13. During the first hour and as a result of the low recharge value (1 × 10 -9 m s -1 ), the water table drops slowly and no erosion occurs. All parameters, except for shear stress, decrease rapidly. An intense storm (recharge of 2.8 × 10 -5 m s -1 ) 545 starts after 1 hr; groundwater seeps out at the cliff base, where the shear stress is highest, forming a gully between the cliff base and the beach. In the ensuing 2 hr, the recharge rate stays the same but the water table height in the sandy area decreases due to the high permeability; as a result the shear stress decreases and erosion stops. At the third hour, the storm stops and there is no recharge; the water table drops slowly and the shear stress and groundwater seepage decrease. The results of the stream power law model with a low intensity rainfall scenario 550 are similar to those of the high intensity rainfall scenario (Fig. 14).

565
The results of the linear diffusion model with high intensity rainfall scenario are shown in Fig. 15. During the first hour of the simulation, the recharge is low (1 × 10 -9 m s -1 ); the water table height, shear stress, diffusivity coefficient and groundwater flux decrease very slowly, whereas the surface water discharge increases. After the 5.9 to 6.5 m), as do all the parameters, reaching their maximum values. A small box canyon develops at the base of the cliff where the shear stress is highest. The intense storm continues for the next 2 hours, and the water table height rises quickly to its maximum value of 7.1 m. The main parameters decrease, and the box canyon enlarges to 16.5 m × 11.3 m, eroding into the cliff. At the end of the storm there is no recharge, the water table drops slowly, the shear stress and the groundwater seepage decrease rapidly, and erosion stops. During the rest of the 575 simulation there is no recharge; the water table continues to drop, the shear stress and the groundwater seepage decrease, and no erosion occurs. The linear diffusion with low intensity rainfall model (Fig. 16)   Finally, the results of the stream power law model without groundwater seepage show a set of closely-spaced gullies, growing from the base of the cliff upwards. The largest of these are incised into the sand unit (Fig. 17).

Discussion
Box canyons are characteristic landforms along the Canterbury coast (Fig. 1a). They are an important driver of coastal geomorphic change as well as loss of agricultural land. The morphologic attributes of the box canyons are similar to those associated with groundwater seepage elsewhere (Table 1).

615
The Canterbury box canyons initiate and evolve via two types of groundwater-related processes. The first process is seepage erosion of sand, which leads to the formation of alcoves and tunnels (Fig. 4). This lowers the overall factor of safety of the slope and is a precursor to the second process, which is slope failure. Box canyons primarily evolve by retrogressive slope failure (Figs. 5,7), which results in the elongation of the box canyons and, to a lesser extent, in widening and branching along the canyon walls. The shear stress of the water seeping out of the canyon 620 head and flowing along the base of the canyon is high enough to entrain the failed sandy material (shear stress of 7.5 × 10 -5 in Fig. 15a; Petit et al. (2015)), which explains why the floors of box canyons are primarily covered by gravel. Wave erosion removes the failed material at the canyon mouths and the base of the cliff. Isotropic scaling of length with width ( Fig. 3c) suggests that canyon planform shape is generally geometrically similar at consecutive stages of evolution.

625
The field observations of box canyon morphology and development are most similar to the LEM results for the high intensity rainfall scenario of the linear diffusive erosion model (Fig. 16). There are two reasons for this: first, as observed in Culling (1963) and Barnhardt et al. (2019), the retrogressive slope failure is best simulated by a diffusive linear model; secondly, our linear diffusive erosion model links the diffusion co-efficient with the shear 630 stress generated by the groundwater. The evolution of the different models and scenarios appears to be influenced by the initial topography -with the cliff preventing the upslope development of erosive features in the stream power law mode, regardless of the recharge values -and the water table height, which depends on permeability and recharge . In the stream power law mode, erosion is linked to surface water, which depends exclusively on groundwater seepage. Since the water table can only intersect the beach or the cliff surface, surface 635 erosion of the plain above the cliff top is not possible . When the stream power law model does not depend on groundwater seepage, the surface erosion does affect the cliff face, but it results in a series of closely spaced and narrow gullies (Fig. 17). The conceptual model, the program configuration and the initial and boundary conditions make the LEM outcomes only applicable to box canyon formation at this site.

640
The location of the box canyons is controlled by two factors. The first factor is the occurrence of sand lenses across a sandy gravel cliff face. This geological framework is conducive to alcove formation, tunnelling and slope failure. The higher permeability of the sand and clean gravel lenses, in comparison to the surrounding sandy gravel (  modelling and inversion, is required to ascertain the sub-surface hydraulic geometry responsible for the alongprofile amplitude variations. The above confirm the importance of spatial variations in hydrogeological properties as a factor controlling the location of a box canyon. This had initially been suggested by Dunne (1990) and has been documented for box canyons in bedrock environments (Laity and Malin, 1985;Newell, 1970). Development

660
of box canyons downslope of permeable conduits may also explain why most of the erosion entails elongation of existing canyons, rather than formation of new ones (Figs. 5, 7). It also agrees with the results of experimental modelling by Berhanu et al. (2012), which suggest that channels grow preferentially at their tip when the groundwater flow is driven by an upstream flow.
Box canyon formation is rapid (daily timescales) and recent (<3 years ago). It is an episodic process that occurs after a threshold is exceeded. This threshold entails a rainfall intensity of >40 mm/day (Figs. 6,7b,11,15), which occurs once every 227 days, on average (Fig. 7b). According to our LEM, up to 100 m 3 of water are estimated to have seeped out of the cliff face to erode 2540 m 3 of material (Fig. 15), which contrasts with the inference by  that 100-1000 times more water than volume of eroded sediment must be discharged in order to create a sapping valley. The erosion rate documented in our study area is up to 30 m per day (Figs. 5e-f), which is the highest rate documented for box canyons so far.
The majority of the box canyons in our study area have shown evidence of erosion in the past 11 years. The OSL data, however, suggest that two largest box canyons have largely been inactive during at least the last 2 ka (recent 675 erosion is only documented in the central section of their mouths, Figs. 5-7). This contrasts with the inference by Schumm and Phillips (1986) that they were formed by spillage of water from swamps behind the cliffs in the 19 th century. Instead we propose that the largest box canyons are relict features that formed as a result of higher groundwater flow, and possibly surface erosion, in the past. The age of sample NZ13A suggests that this may have occurred during the Last Glacial Maximum.

6 Conclusions
Box canyon formation is a prevalent process shaping the Canterbury coast of the South Island of New Zealand.
In this study we have integrated field observations, OSL dating, multi-temporal UAV and satellite data, time-685 domain electromagnetic surveying, and slope stability and landscape evolution modelling to constrain the controlling factors and temporal scales of box canyon formation. Our results indicate that box canyon development in sandy to gravel cliffs is a groundwater-related, episodic process that occurs when rain falls at intensities of >40 mm per day. At the study area, such rainfall events occur at a mean frequency of once every 227 days. Box canyons have been developing, primarily by elongation, in the last 11 years, with the latest episode dating to 3 690 years ago. Canyons form within days and erosion rates can reach values of up to 30 m per day. The two largest box canyons in the study area, however, have not been actively elongating or widening in the last 2 ka. The key processes responsible for canyon development are the formation of alcoves and tunnels in sandy lenses by groundwater seepage erosion, followed by retrogressive slope failure. The latter is a result of undermining and a decrease in shear strength due to excess pore pressure development in the lower part of the slope. The location of 695 the box canyons is controlled by the occurrence of hydraulically-conductive zones, which comprise relict braided river channels and possibly tunnels, and sand lenses exposed across the sandy gravel cliff. Box canyon formation is best represented by a linear diffusive model and geometrical scaling.