Can artesian groundwater and earthquake-induced aquifer leakage exacerbate the manifestation of liquefaction?

Vast quantities of liquefaction ejecta repeatedly inundated properties during the 2010-2011 Canterbury earthquake sequence in New Zealand, resulting in differential ground surface subsidence and significant damage to buildings and urban infrastructure. There are strong spatial correlations between the occurrence of ejected sediment with groundwater pressure in deep aquifers. When geotechnical testing sites are grouped according to liquefaction vulnerability indices (to control variance relating to shaking strength, water table depth, and soil profile strength), places where ‘minor’ and ‘moderate-severe’ liquefaction occurred during the Mw6.2 Christchurch earthquake had distinctly higher aquifer pressure than sites where liquefaction was not observed. Together with observations of earthquake-induced pressure changes and inferred transfer of groundwater from deep aquifers to shallower levels, an interpretation is that leakage from aquifers with artesian (above ground) pressure provided an additional driving mechanism for surface manifestation of water and sediment. It is surmised that above-ground aquifer pressures further promoted suffusion and piping along fractures, flow-pathways and liquefied horizons. The Mw6.2 Christchurch earthquake is presented as an example where leakage of artesian groundwater likely contributed to the near-surface liquefaction-induced ground damage. The process can result in underprediction of liquefaction vulnerability so needs to be considered when evaluating potential for earthquake-induced liquefaction and ground damage wherever groundwater is confined.


Motivation
Liquefaction has been observed during almost all moderate-to largeearthquakes in New Zealand (Saunders and Berryman, 2012 and refs therein), so the potential for liquefaction in Christchurch City during a large earthquake had been anticipated due to the co-occurrence of a near-surface water table and low-to medium-density non-plastic soils (Brown and Weeber, 1992;Brown et al., 1995). Shaking from a protracted series of seismic events during 2010-2011 triggered widespread and repeated ejection of sand, silt and malodorous sediment-rich water in eastern Christchurch Orense et al., 2011;Quigley et al., 2013;van Ballegooy et al., 2014b;Lees et al., 2015). Groundwater discharged from fissures and eruptive centers, some growing sufficiently large to engulf vehicles and undermine buildings, and flooded properties and streets (Supplementary Materials Movie S1). As water recessed, drained away and evaporated, sandy-silt dried into semi-solid congealed deposits, covering streets and properties, that in places exceeded 0.5 m thickness. The liquefaction-induced ground damage involved a combination of subsidence, differential settlement, lateral spreading, and sediment invasion. It rendered about 15,000 homes beyond economic repair, impaired around 50,000 of the 140,000 residential properties, and damaged approximately 50% of the horizontal infrastructure (roads, electricity, conveyances for waste water and fresh water) (Rogers et al., 2014;van Ballegooy et al., 2014b). Economic losses exceeded US$21 billion, or 15% of New Zealand's GDP, making it one of the world's largest and most-complicated natural disaster insurance events (King et al., 2014;Rogers et al., 2014). It also had major impact on societal health and wellbeing (Potter et al., 2015). With poor ground performance a key factor, the earthquake sequence was one of the most pervasive and extensive urban liquefaction events ever experienced.
This study investigates whether the groundwater in aquifers beneath Christchurch and deep pressure increases caused by shaking could have played a role in exacerbating the liquefaction-related damage. Building on published evidence that the Canterbury earthquake sequence changed sediment permeabilities (Rutter et al., 2016;Weaver et al., 2019) and locally induced upwards flow of groundwater (Cox et al., 2012;Gulley et al., 2013;Dudley Ward, 2015), we explore whether the extensive ejection of groundwater and sediment in the city might be linked to an underlying aquifer system with substantial confined artesian (above ground) groundwater pressure. Collating hydrological, geological, shaking and geotechnical data, we defined the hydrological pressure regime during the Mw7.1 Darfield (4 September 2010, 04:35 NZST) and Mw6.2 Christchurch (22 February 2011, 12:51 NZDT) earthquakes, then looked at relationships between hydrological parameters, mapped occurrence of ejected sediment, and geotechnical parameters used to evaluate site liquefaction vulnerability.
Analysis has been carried out in two parts at different scales. Firstly, spatial correlations are assessed across a wider study area straddling the transition from inland semiconfined to coastal confined aquifers, to provide contextual setting over a wide variety of hydrogeologic conditions and groundwater pressures. Spatial probability distributions indicate strong relationships between the occurrence of ejected sediment and deep aquifer pressure. However, as there is always potential for spatial correlations to be coincidental associations rather than due to some direct causal relationship, the second part of our examination focused on a more-restricted urban area where geotechnical testing data are available. Albeit more biased towards places where damage occurred and within a more-restricted range of hydrogeological conditions, liquefaction vulnerability indices from geotechnical tests provide local proxies that account for combined variations in shaking, soil profile strength, and water table depth. Geotechnical sites were grouped and examined according to their expected liquefaction vulnerability and/or soil behavior-type, irrespective of their spatial location, to explore any systematic variation with deep aquifer pressures below the testing site. Although we consider one setting, the paper collates evidence that aquifer systems can be perturbed and damaged by earthquakes, with potential for leakage to exacerbate liquefaction ejection and surface manifestation.

Hydrogeological setting
The Canterbury Plains are a coalesced alluvial fan built by Quaternary outwash of gravel and sand eroded from the Southern Alps ( Fig. 1) (Brown and Weeber, 1992;Brown et al., 1995). They host an aquifer system that supports nationally important agriculture and is the principal source of drinking-water for about 350,000 residents (Brown, 2001). Groundwater, sourced from river seepage and rainfall infiltration, flows seaward down the topographic gradient from mountains towards the coast (Fig. 1B). Subsidence through the Quaternary period, at long-term rates of 0.25-0.44 mm/a (Begg et al., 2015), and relative sea level fluctuations during glacial-interglacial cycles, has resulted in a coastal zone where gravels are inter-bedded with wedge-shaped estuarine and shallow marine sequences. Channelized gravels re-worked and sorted by alluvial processes form particularly good aquifers whereas finer-grained marine/estuarine silt, silty-sand, clay, and peat layers act as aquitards that can be found up to 15 km inland from the present-day shoreline (Brown et al., 1995). Low-permeability fine-grained sediments restrict vertical flow and start to confine groundwater at thicknesses of >1-3 m (Talbot et al., 1986;Weeber, 2008). The shallowest confining estuarine/marine units (Christchurch and Upper Springston formations) reach around 40 m thickness beneath the present-day coastline (Fig. 2).
Christchurch City lies at the coastal edge of Canterbury Plains, straddling the transition where groundwater of inland semi-confined and unconfined aquifers becomes progressively confined beneath marine/estuarine sediments (Figs. 1, 2). West of the city, the depth to groundwater exceeds 10 m and a downwards component of flow (negative vertical head gradient) recharges groundwater (Brown and Weeber, 1992;Weeber, 2008; van Ballegooy et al., 2014a). The western suburbs of Christchurch are constructed on paleo-channel and overbank deposits left by rivers that once drained the Southern Alps, with depth to groundwater typically 5-10 m. Springs emerge where groundwater flow becomes impeded by fine-grained sediment, feeding the Avon/Ōtakaro and Heathcote Rivers that flow through the city. Beneath eastern Christchurch, a wedge of marine/estuarine sediments (Christchurch and Upper Springston formations) hosts a shallow unconfined water table at <2 m depth that is hydrologically connected to the rivers (Fig. 2C, D), overlying at least four confined deep gravel aquifers (being the Riccarton 20-40 m; Linwood 50-80 m; Burwood 100-110 m; and Wainoni 130-150 m gravels). The hydraulic head of deep aquifers beneath eastern Christchurch reaches 5 m or more above-ground and vertical head gradients are upwards (positive) with geopressure consistently 5%, locally exceeding 10%, above hydrostatic pressure (Fig. 2E). A small component of vertical flow seeps through the aquifer system, although exact volumes have been difficult to quantify (Lough and Williams, 2009). A series of new springs emerged through confining sediments following the Mw7.1 Darfield (4 Sept 2010) earthquake, maintaining sustained flows that became locally problematic (Cox et al., 2012). The earthquakes caused changes in aquifer transmissivity that induced sustained groundwater rises near Greendale fault, but the hydraulic head in confined aquifers beneath Christchurch was lowered by around 1 m (Rutter et al., 2016).

Data and methods
This study examined the two major events of the Canterbury earthquake sequence that induced widespread liquefaction-related land damage: the Mw7.1 Darfield (4 September 2010, 04:35 NZST) earthquake and the Mw6.2 Christchurch (22 February 2011, 12:51 NZDT) earthquake. Aftershocks of Mw5.5 (13:05) and Mw5.6 (13:05) closely followed the Mw6.2 earthquake on the 22 February, so observations of liquefaction and groundwater pressure changes could not always be distinguished into separate shaking episodes during this earthquake event (Quigley et al., 2013(Quigley et al., , 2016see below). Several other earthquakes caused liquefaction, but progressive damage to the groundwater monitoring network, or close timing of foreshocks and aftershocks, makes the groundwater behavior during these latter events even more difficult to constrain. We collated groundwater, geological, shaking and geotechnical data within a 54 × 40 km study area (172.25 • E, 43.25 • S to 172.75 • E, 43.75 • S) that surrounds and includes Christchurch City, using some external data to extend interpolations right to the study area boundary. We then focus on a more-limited 25 × 22 km urban region, where liquefaction was most prevalent and many geotechnical tests have been completed (Fig. 1A).

Groundwater monitoring data
Groundwater monitoring by Environment Canterbury (ECAN), Christchurch City Council (CCC), and the Earthquake Commission (EQC) provides a comprehensive public dataset spanning many decades that includes the earthquake sequence ( Fig. 2A,B,C). ECAN has been mostly focused on understanding groundwater in deeper aquifers. There Christchurch urban area, inland semi-confined and unconfined aquifers (yellow), coastal zone of confined aquifers (green), piezometric groundwater contours (5 m intervals above sea level), and epicenters of the Mw7.1 Darfield (4 Sept 2010) and the Mw6.2 Christchurch (22 Feb 2011) earthquakes. (B) Schematic block diagram of the aquifer system and groundwater flow regime (adapted from Weeber, 2008). (C) Potentiometric levels relative to ground surface, highlighting the west-east difference in head between the shallow water table compared with deeper aquifer groundwater. Note that colour figures are available in the electronic version of this paper. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2.
Hydrogeology and geopressure. (A) Map of wells with weekly or monthly data used to develop grid models of static pre-earthquake piezometric groundwater levels, coded according to different geological formations (Brown and Weeber, 1992;Brown et al., 1995). The legend orders formations from shallowest to deepest. (B) Distribution of ECAN monitoring wells measuring aquifer pressure at 15-min intervals used to define dynamic head changes during the earthquakes. (C) Grid model and contours of the water table depth relative to ground surface in February 2011. Small crosses are wells and piezometers that constrain the model between rivers. (D) Isopach model of shallow aquitard thickness (refined after Weeber, 2008). Grey areas correspond to 0 m depth to gravel/confining sediment thickness. (E) Graph of hydraulic head relative to local ground compared against depth of wells, colored according to geological formation (as per Fig. 2A). (F) Map of vertical head gradients modelled during the Mw6.2 22 Feb 2011 earthquake, based on the difference between the total coseismic aquifer piezometric level and the local elevation of the water table, divided by the depth between the water table and the bottom of the well screen. Positive head gradients with upwards flow are orange-red, and downwards negative flow are light-dark blue. The green contour outlines where total coseismic head in deep aquifers are 0 m relative to ground. Maps have black lines showing the Greendale (easternmost end -solid) and Port Hills (concealed -dashed) faults that ruptured during the earthquakes, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) are ~ 170 wells across Canterbury where weekly-monthly observations are made and ~ 76 other wells hosting automated transducers recording at 15-min intervals. There were ~ 40 wells with transducers recording in the study area during the Mw7.1 Darfield earthquake (Fig. 2B), but the deep monitoring network suffered attrition during the earthquake sequence. By way of contrast, CCC and EQC monitoring characterizes the phreatic water table within the Christchurch urban area where >960 shallow piezometers (depths <10 m) had monthly manual observations from 2010 to 2013 as the network was developed for earthquake recovery (Fig. 2C). ECAN also holds sporadic data from one-off surveys, consent measurements, and irregular observations in their database; these data are useful for constraining interpolation between automated monitoring sites. River levels have been monitored at 12 sites within the study area and were used, together with surveyed stream channel profiles, to locally constrain the water table elevation (cf. van Ballegooy et al., 2014a).
Rather than converting groundwater level measurements to an absolute elevation, we deliberately analyzed and report all groundwater relative to local ground level, which avoids compounding well-head surveying elevation uncertainties into calculations. Groundwater pressure (absolute relative to ground) or pressure changes presented as hydraulic head, or changes in head, in units of meters. Calculations assumed a freshwater composition at 25 • C and density 0.997 g/mL. Well data were first imported into a 3D geological model of Christchurch region in Leapfrog® (Begg et al., 2015) where previous classifications by ECAN were reviewed on the basis of (i) the geological formation in which they were screened ( Fig. 2A) and (ii) aquifer-type based on well depth, hydraulic head relative to ground, and data from nearby wells (Fig. 2E). Isopach contours of confining sediment thickness, inferred from manual interpretation of the depth at which gravel was first encountered in the ECAN well logs (Weeber, 2008), were also evaluated in the 3D geological model (Begg et al., 2015). These were slightly refined, principally because of new drilling and the ability to look for anomalous points in the 3D model, then interpolated to produce an aquitard thickness map (Fig. 2D) using the 'Topo to Raster' function in ArcGIS (a discretized thin-plate spline technique that is specifically formulated to work with co-dependent data inputs).
Specific details of hydrological responses of individual wells to the earthquakes, which include a variety of changes in water level and tidal sensitivity, have been published elsewhere (Cox et al., 2012;Gulley et al., 2013;Rutter et al., 2016;Weaver et al., 2019). As well as interpolating maps that highlight and summarize the dynamic pressure changes, we generated movies depicting the rise and fall of groundwater in monitoring wells and subsequent recession to new post-earthquake levels (Movies S2, S3 Supplementary Material). They show 15-min monitoring data at each well as a function of distance from the earthquake epicenter, all with immediate increases following each earthquake, then establishing new equilibria over subsequent hours to days, at reduced pressures in deep wells or elevated pressures in shallow wells. Previous work concludes that the well responses provide evidence for earthquake-induced transfer of water upwards from deep to shallow levels (Cox et al., 2012;Gulley et al., 2013;Dudley Ward, 2015). Modelling of groundwater level changes suggests that upwards flow involved both an increase (12-19%) in the relative conductance of a confining aquitard and decreases in sediment storativity (by 5% in a shallow unconfined aquifer and 2% in deep confined aquifers) (Dudley Ward, 2015). The earthquakes also altered the tidal influence on some confined aquifers (Weaver et al., 2019), and the shallow unconfined water table remained anomalous for days to weeks after each earthquake (van Ballegooy et al., 2014a;Quigley et al., 2016). For this study, we calculated 'static' pre-earthquake piezometric levels for all wells based on the mean recorded values during the month immediately prior to the Mw7.1 Darfield (4 Sept 2010) and Mw6.2 Christchurch (22 Feb 2011) earthquakes, and the 'dynamic' pressure changes induced by the shaking at monitoring wells with 15-min data records following the method of Gulley et al. (2013).

Interpolation of static and dynamic groundwater pressures
Groundwater potentiometric surfaces were derived by interpolating data between wells into grids with 50 × 50 m cells using the natural neighbor method in ArcGIS (Sibson, 1981). This is a local polynomial that produces a smooth surface that ensures interpolated points lie within the range of original observations. 486 wells and piezometers screened in the shallow unconfined water table (mostly Christchurch and Upper Springston formations) were distinguished from 630 wells in 'deep' confined and semi-confined aquifers (the Lower Springston and Burnham Formation, or Riccarton, Linwood, Burwood, and Wainoni gravels). Initial attempts to develop individual surfaces for each aquifer were impaired by the small number and/or uneven distribution of wells, so we simplified to a combined 'deep' aquifer model and distinct phreatic, shallow 'water table' surface. By treating all confined aquifers (Lower Springston, Riccarton, Bromley, Linwood, Burwood, and Wainoni Formations) together as 'deep' we made a trade-off, whereby the final surfaces contain some local anomalies in absolute values, but the spatial resolution and map interpolation is greatly improved for calculated parameters -such as the vertical head gradient (Fig. 2F). Interpolation of the unconfined water table (as a depth to groundwater in the Christchurch, Upper Springston, and Riccarton Formations) was also constrained using river gauging data with channel surveys and the position of the mid-tide sea level at the coast (following van Ballegooy et al., 2014a).
'Static' pre-earthquake hydraulic heads in confined-and semiconfined aquifers were interpolated from 'deep' wells using the mean values for the month preceding the earthquakes (Fig. 3A,D). 630 wells in the study area have one or more water-level observations during Aug-Sep 2010. There were far fewer with continuous records during the earthquakes: monitoring at 15-min intervals occurred at 42 'deep' sites during the 4 Sept 2010 Mw7.1 earthquake, but only 25 sites during the Feb 2011 Mw6.2 due to network damage and attrition. 'Dynamic' head changes were derived for the earthquake-induced groundwater level changes from observations immediately before, and peak pressure after, each earthquake according to the method of Gulley et al. (2013) (Fig. 3B,C). There were 5 (of 25) wells where groundwater pressure changes on 22 Feb 2011 had an additional contribution of pressure from the Mw5.6 aftershock at 14:50, before the response to the Mw6.2 earthquake at 12:51 had fully dissipated. As aftershock contributions are not easily distinguished in these data, the entire pressure response was attributed to the Mw6.2 earthquake. Total 'coseismic' hydraulic head for 'deep' confined and semi-confined aquifers were then modelled by adding interpolated grids of dynamic head change to the static preearthquake piezometric levels ( Fig. 3E,F).
Dimensionless vertical groundwater head gradients were also interpolated from site calculations. The difference between the piezometric level in each well and the local depth to water table (ΔH in m) was divided by the vertical distance between the water table and the lowermost elevation of the monitoring well screen (Δz in m). We present and analyze only the coseismic vertical head gradients modelled for the Mw7.1 and Mw6.2 earthquakes (e.g. Fig. 2F), as the patterns of static, pre-earthquake, upward, and downward gradients are similar, but less directly relevant to conditions during the earthquake. We also modelled depth to the unconfined shallow water table (Fig. 2C) and an isopach map of shallow confining sediment thickness (Fig. 2D).

Land damage mapping
Qualitative surveys of land damage were carried out immediately after major events in the Canterbury 2010-2011 earthquake sequence. In the urban areas of Christchurch and the northern township of Kaiapoi, detailed ground-based inspections recorded land and dwelling foundation damage with a measure of severity, using five or six classes ranging from no observed sediment ejection or cracking to large quantities of ejected material and severe lateral spreading (Tonkin and Taylor Ltd, 2013; van Ballegooy et al., 2014b). Observations were recorded onto land-parcel polygons averaging 12,500 m 2 , with a few >100,000 m 2 . Non-residential land parcels were not assessed by these ground surveys, so areas of sediment ejected in streets, parks, and some business areas were mapped separately from high resolution aerial photographs collected shortly after each earthquake. The quantity of ejected sediment classified as 'none', 'minor', or 'moderate-severe' (Canterbury Geotechnical Database, 2015). For the rural area surrounding Christchurch, mapping from high-resolution aerial photographs and satellite imagery (WorldView 1 & 2) identified places where there was manifestation of sediment, flooding or sediment-laden flood water, together with a measure of observational confidence (Townsend et al., 2016).
For this study, we combined the available urban and rural land damage datasets, utilizing only those remote observations of ejected sediment that had been assigned as 'certain' or 'probable', ignoring observations that were 'possible' or 'uncertain'. We then classified the ejected sediment into simplified 'not observed' or 'occurs' classes, and where present in urban areas with sub-classes of 'minor' or 'moderatesevere' (volume or severity) (Fig. S1). Land-parcel polygon maps were transformed into continuous 50 × 50 m raster grids to match interpolated hydrological grids, covering the 2160 km 2 area of interest with a maximum 864,000 cells. There was surface manifestation of sediment over 12% of the Christchurch City's developed land following Mw7.1 Darfield earthquake (4 Sept 2010) earthquake, and 24% following the closer Mw6.2 Christchurch event (22 Feb 2011), which corresponds to ~ 4% and ~ 9% of our wider study area, respectively.
There are several important caveats with regard to the land damage mapping. Firstly, the maps only identify features visible at the time of the inspection/photography and cannot delineate where liquefaction processes occurred in the subsurface, or places where clean-up or changes might have occurred prior to inspection. For this reason, we deliberately adopt the terminology 'not observed' over 'not present' to describe the absence of ejected sediment and distinguish 'ejected sediment' as observed material that is distinct from the process of 'liquefaction'. We chose not to include areas of flooding observed in aerial photographs and satellite imagery in our analyses, because there had been some rain, evaporation, infiltration and lateral flow of ejected water between the time of the earthquake and the capture of some images. Land damage and its mapping may in part be cumulative, in that observations following the Mw6.2 earthquake on 22 Feb 2011 potentially include sediment ejected by the Mw7.1 earthquake on 4 Sept 2010 and/or immediate aftershocks, although liquefaction-induced ground damage from the Mw6.2 earthquake was greater and more widely distributed. Following earthquakes in June and Dec 2011 it was not always possible to distinguish newly ejected sediment and fresh damage from that caused previously during larger events or closely spaced earthquakes; another reason for focusing this study on the two main (Mw7.1 and Mw6.2) earthquakes.

Geotechnical data
The Canterbury Geotechnical Database (CGD) was developed during the recovery of Christchurch as a collaborative sharing model for public and private data that is internationally unique (Scott et al., 2015). It has since grown into the New Zealand Geotechnical Database (NZGD) where Christchurch data are available (New Zealand Geotechnical Database, 2020). With extensive geotechnical investigations commissioned by the New Zealand Earthquake Commission (EQC) and private insurers in areas affected by consequential liquefaction, there are data from over 20,000 cone penetration (CPT) tests and 1000 piezometers in Christchurch that enable close examination of ground strength and performance. Although there is a large quantity of geotechnical information available, it is significantly biased towards sites where liquefaction manifestation occurred in the Christchurch urban area where assessments were needed for recovery or insurance purposes. There are fewer test data available in places where there was no liquefaction.
Geotechnical testing was centered on the Christchurch urban area and immediate surrounds, covering a 25 × 22 km area (but totaling 325 km 2 ) within the wider 54 × 40 km area considered by this study, almost entirely within the coastal zone of confined aquifers (Fig. 1A). Testing is shallow, typically to depths of 10 to 15 m and rarely exceeding 20 m, following government guidelines (MBIE, 2012(MBIE, , 2015. Although most testing was completed after earthquakes, there does not appear to be any significant change in measurable ground strength as the earthquake sequence evolved (Orense et al., 2012;Lees et al., 2015;Russell et al., 2015 Appendix E). The limited distribution of testing relative to changes in aquifer geology, and the bias towards testing sites where liquefaction manifestation occurred, was a key reason for presenting this study in two parts at two scales.
Engineering assessments of liquefaction vulnerability commonly involve semi-empirical or advanced numerical methods to derive parameters that define triggering of the liquefaction process (e.g. Iwasaki et al., 1978;Seed and Idris, 1971;National Research Council, 1985;Robertson and Wride, 1998;Youd et al., 2001). Most approaches derive a cyclic stress ratio (CSR), representing the seismic demand on a soil layer caused by an expected earthquake, and compare it against a cyclic resistance ratio (CRR), representing the capacity of the layer to resist liquefaction. The transformation from a solid to liquefied state is due to the increase in pore water pressure that accompanies cyclic deformation of non-cohesive sediment. A total loss of stiffness and strength occurs when pore water pressure exceeds the effective overburden stress of the material (Seed and Lee, 1966;Terzaghi et al., 1996). We utilized geotechnical parameters from CPT testing calculated by van Ballegooy et al. (2015b), removing any soundings where there was greater than 3 m of predrilling, and any tests that did not exceed depths of 5 m. Although published parameters from Christchurch also include settlement deformation (S) (MBIE, 2015) and the Liquefaction Severity Number (LSN) (Tonkin and Taylor Ltd, 2013; van Ballegooy et al., 2014b), for simplicity we report only the Liquefaction Potential Index (LPI) (Iwasaki et al., 1978). We also utilize an average Soil Behavior Index (Ic10) determined over the top 10 m of each CPT profile.
The LPI calculations from van Ballegooy et al. (2014a) used the IB-2014 liquefaction triggering procedure (Idriss and Boulanger, 2014), peak ground acceleration (PGA) models of Bradley and Hughes (2012), and event-specific interpolated depths to the water table (Tonkin and Taylor Ltd, 2013). We used 12,005 LPI site values determined for the Mw7.1 earthquake and 11,603 for the Mw6.2 earthquake, which reflects differing spatial extent of the input models. There were marked differences in the shaking experienced across the study area during the Mw7.1 and Mw6.2 earthquakes (Bradley, 2014;Bradley et al., 2014). These differences are highlighted in maps of PGA that show attenuation from southwest-northeast and lower (~0.15 to 0.25 g) accelerations in the Christchurch urban area during the Mw7.1, and attenuation from southeast-northeast and higher (~0.25 to >0.5 g) accelerations in the urban area during the Mw6.2 (see Fig. S1C,D). Correspondingly, LPI values calculated for the Mw7.1 earthquake were generally lower than those for the Mw6.2 event.
Several studies have shown dispersion between predicted and observed liquefaction in Christchurch, which are typically attributed to uncertainties in ground motion and site characteristics, anthropogenic modifications and infrastructure, and/or use of simplified models (e.g. Maurer et al., 2014;van Ballegooy et al., 2015b;Cubrinovski et al., 2019;Borella et al., 2020). The degree of dispersion is only marginally affected by parameter choice (see van Ballegooy et al., 2015b). We utilized locally derived threshold values for classifying LPI into categories of low (LPI ≤ 5), moderate (5 < LPI ≤ 15) and high (LPI >15) vulnerability (after Maurer et al., 2014;Wotherspoon et al., 2015 and refs. therein). Where the Ic10 provides an indication of soil behavior type, the LPI values are a more general proxy that combine all variations in seismic load, water table depth and shallow soil strength profile. We utilize geotechnical data in two ways: either (i) as individual site values, with observations and data attributed to that site from other datasets, under the assumption that there is independence between sites; or (ii) as interpolated 50 × 50 m grids, for correlating against land damage mapping and interpolated hydrological data in assessments of spatial correlation (see Supplementary Materials). By interpolating site data, we make an implicit assumption that there is some degree of lateral codependence in the liquefaction vulnerability between nearby sites, that might be caused by common subsurface geology, hydrology, or peak ground acceleration. We used the natural neighbor method (Sibson, 1981) to interpolate geotechnical data. Where two or more geotechnical tests were < 50 m apart, the cell was assigned the 'worstcase' highest calculated LPI value.

Methods
Groundwater parameters show many spatial similarities with maps of ejected sediment in Christchurch (Figs. 4A, S1) which led us to explore a general hypothesis "that groundwater pressure in confined aquifers exacerbated the liquefaction hazard". To further quantify the spatial correlations between the aquifers and surficial damage we derived histograms and spatial probability distributions. Each variable was assessed in bins of equal interval, then the number of 50 × 50 m cells of each ejected sediment category ('occurs' vs 'not observed'), or severity subcategory ('minor' or 'moderate-severe'), within each bin range was counted and converted to histograms (Fig. 4B). Spatial probability distributions were calculated by determining the proportion of each category or subcategory within each separate bin independently of other bins (Fig. 4C). In essence each bin range is treated as a simple binomial (or multinomial) experiment, where the number of trials is equal to the total number of cells in the bin range, with a probability derived from the number of occurs 'successes' (or 'minor' vs 'moderatesevere' successes). On average, liquefaction was manifest over ~ 4% of the 54 × 40 km study area during the Mw7.1 Darfield earthquake and ~ 9% during the Mw6.2 Christchurch earthquake. The bins have a lower total number of cells and larger uncertainty in occurrence 'success' probability towards the maximum and minimum tail of each variable. Uncertainties across the probability distribution were approximated using standard errors, which are inversely proportional to the square root of the sample size. A set of cumulative curves were then calculated for each ejected sediment category using the number of 50 × 50 m cells within each bin interval as a percentage within the total population of each sediment category, treating each category independently.
The probability distributions were used to derive observations, refine hypotheses and develop mechanistic explanations. We looked for positive or negative correlations between ejected sediment and the different groundwater, geological, shaking and geotechnical parameters, particularly where counts suggested high probabilities associated with a particular variable (e.g. >50%). Results for the Mw6.2 Christchurch had many similarities to those from the Mw7.1 Darfield earthquake, but since surface manifestation of sediment was more pronounced during the Mw6.2 (Fig. S1) the trends from the smaller, more proximal and devastating earthquake have lower uncertainty and appear more pronounced. The figures selected for the paper are focused towards the Mw6.2 event, but equivalent plots for the Mw7.1 are provided in supplementary material (e.g. Fig. S2).
To examine and test hypotheses, we apply a number of empirical tests to both site and grid datasets. Population statistics are plotted with median lines, boxes spanning to 25th and 75th percentiles, triangles to 10th and 90th percentiles, and 'whiskers' spanning the full data range. Box widths are scaled proportional to the total sample in each plot or plot panel, normalized such that the sum of the widths in each plot partition is constant. Two sample Kolmogorov-Smirnov (KS) tests were also carried out using p-values to test degrees to which populations are significantly different. But instead of collapsing data down to tables with single p-value variables and/or represented by linear regression, we present a series of box and whisker plots as a means to visually assess population differences.

Spatial probability distributions for ejected sediment versus groundwater
Occurrence of ejected sediment compared with the depth to water table, binned at 0.2 m intervals (Fig. 5A), showed a strong spatial relationship during the Mw6.2 earthquake: ejected sediment occurred in 24% of the study area where the water table was at surface (<0.2 m); where the water table was deeper the proportion of observed sediment ejection decreased steadily; from 4 m or deeper there were very few occurrences of ejected sediment (<1%). Surface manifestation has spatial correlation with water table depth, but the severity of manifestation may not (see also van Ballegooy et al., 2014bvan Ballegooy et al., , 2015a. Results during the more-distal Mw7.1 earthquake had many similarities (Fig. S2A), but ground accelerations were lower and there were fewer places where the water table was shallower than 3 m depth during this spring-season earthquake. That 'occurrence' could be related to groundwater depth was an expected result (Ishihara, 1985;van Ballegooy et al., 2015b).
Dynamic head changes of up to 5.7 m (~ 57 kPa) were recorded in deep aquifers during the Mw6.2 Christchurch earthquake, but were locally much less-pronounced in the study area (up to 2.1 m, ~ 21 kPa) during the more-distal Mw7.1 Darfield earthquake (Fig. 3B,E; Movies S2, S3). Aquifer pressure changes are broadly correlated with well depth and earthquake magnitude/distance (Cox et al., 2012;Gulley et al., 2013), being largest beneath the city in wells screened in the deepest, confined coastal-zone aquifers. Artesian wells subsequently recorded post-earthquake recovery to lower pressures (negative offsets), then remained at ~1 m lower pressure for a few years (Rutter et al., 2016), indicating that the earthquakes induced a loss of groundwater from deep aquifers to shallower levels (Cox et al., 2012;Gulley et al., 2013;Dudley Ward, 2015). Mapped occurrence of ejected sediment has a strong spatial correlation with the dynamic pressure change during the Mw6.2 earthquake, such that where head increases exceeded 3 m there was over 50% surface manifestation of sand and silt (Fig. 5B). Ejected sediment was only observed where some form of positive head change (>0 m) was experienced. During the more-distal Mw7.1 Darfield earthquake, the dynamic head change was smaller and less variable across the study area, and did not produce as clear a correlation with the mapped distribution of ejected sediment (Fig. S2B). Some moderate-severe sediment ejection and ground damage did occur in Kaiapoi and eastern Christchurch during the Mw7.1 earthquake in the absence of any dynamic change in head.
Relationships between ejected sediment and the total coseismic head in deep aquifers, binned into 1 m intervals relative to ground surface (Figs. 4,5C; Fig. S2C), show that where the total coseismic head remained lower than − 2 m (i.e. > 2 m below ground) during the Mw6.2 Christchurch earthquake there were very few occurrences of sediment manifestation (<1%); as pressures became positive (artesian) occurrences of sediment rose sharply; where coseismic head exceeded 6 m, over 50% of the area had surface manifestation, rising to over 75% occurrence at heads of 9 m.
Of the hydrogeological parameters tested in this study, the spatial probability distribution for ejected sediment against total coseismic head during the Mw6.2 earthquake (Fig. 4,5C) was the strongest spatial relationship observed. Its apparent ability to predict the manifestation of liquefaction appears comparable to LPI (see below and Fig. S3B) in this instance. There are many similarities between the spatial probability distributions for the more-distal Mw7.1 earthquake (Fig. S2) and equivalents for the Mw6.2 (Fig. 5), but the total area of ejected sediment during the Mw7.1 was lower (4% in the study area), peaked at 41% (where total coseismic head was 6-7 m), and probability distribution cumulative curves for severity subclasses are not as clearly distinguished.
Relationships between ejected sediment occurrence and the vertical head gradient or confining sediment thickness hint at possible contribution of the confined aquifer system to processes causing manifestation of liquefaction. Ejected sediment occurrences nearly all coincided with places where the coseismic vertical head gradient was positive (Fig. 5D). An interpretation is that a positive gradient, and/or upward flow of groundwater, promotes movement of sand and silt-sized sediment towards the surface. Where vertical head gradients are positive, and hydraulic connections exist between near-surface liquefied sediment and underlying groundwater at higher pressure, shallow porepressures already affected by cyclic deformation can potentially be further increased (Wang, 2007;Wang and Manga, 2010;Cubrinovski et al., 2019). Spatial probability distributions of ejected sediment occurrence versus thickness of confining sediment are positively correlated with isopach thickness (Fig. 5E): albeit noisy, rising gently to over 20% of areas where the aquitard is 20 m thick, and more steeply to reach ~ 50% where it is over 30 m thick. The marine-estuarine silty-sands and sandy-silts may be more prone to sediment ejection in areas where they are thicker, less gravelly, and younger. But it is difficult to isolate whether the spatial correlation is causative or merely coincidental, as the increasing aquitard thickness covaries with two other factors controlling liquefaction: the southwest to northeast shallowing of the unconfined water table (Fig. 2C) and the northwest-southeast increase in shaking (peak ground acceleration) with proximity to the Mw6.2 earthquake epicenter (Fig. S1D). Total coseismic hydraulic head also increases as the aquitard thickness increases from northwest to southeast. Since it is plausible that spatial correlations are coincidental rather than mechanistic, we explore the possibility further in the next section.

Nonspatial analysis of explanatory variables
There are several spatial variations across the study area that could produce coincidental associations. Water table depth, confining sediment thickness, sediment facies and grainsize (from non-marine to marine), and deep aquifer pressure all show general change from westeast. There were also gradients in seismic loading (peak ground acceleration in particular) during the earthquakes. To examine relationships and causative explanations further, we now focus on the urban region, where liquefaction was more prevalent and there are many geotechnical tests. Although the groundwater parameters have more restricted ranges in the urban area, and tests are biased towards places with liquefaction damage, calculated LPI values (van Ballegooy et al., 2015b) at ~12,000 (caption on next column) Fig. 5. Ejected sediment vs. hydrogeology, Mw6.2 Feb 2011. Spatial probability distributions and cumulative percent curves for ejected sediment occurrence (including minor and moderate-severe subclasses in the urban area) and 'not observed' compared against hydrogeological parameters. The spatial probability distributions show the % of cells of each category within each bin, whereas cumulative distributions are calculated from the number of cells normalized within the total population of each category independently. (A) Water table depth (below ground) binned at 0.2 m intervals. (B) Dynamic aquifer head change induced by the earthquake, binned at 0.1 m intervals (as per Fig. 3E). (C) Total coseismic hydraulic head (as per Fig. 3F). (D) Vertical head gradient in 1% intervals relative to normal hydrostatic pressure. (E) Thickness of confining sediments (Christchurch and Upper Springston Formation) in 1 m intervals. sites incorporate measures of ground strength, earthquake shaking intensity and soil saturation (inferred by the shallow unconfined water table position) to provide a prediction of the likelihood of liquefaction (see also Fig. S3). We start by evaluating the null hypothesis 1 "that observed occurrence or severity of ejected sediment is unrelated to liquefaction parameters when controlled for aquifer pressure" through ordinal multinomial responses developed between site liquefaction indices, ground damage mapping, and the groundwater surfaces.
First, we analyzed the liquefaction potential index (LPI) against the categories of observed ejected sediment ('not observed', 'minor', 'moderate-severe') for each earthquake, with sites categorized by the total coseismic hydraulic head in the underlying aquifers (< 0 m below ground, 0 ≤ head ≤ 6 m, head > 6 m) (Fig. 6A,B). There are no clear trends between liquefaction occurrence and LPI values during the Mw7.1, or differences from one aquifer pressure regime to another (Fig. 6A). Within each aquifer pressure category during the Mw6.2, there was a tendency towards higher median LPI values at places where ejected sediment was observed (either minor or moderate-severe) compared with not observed (Fig. 6B). Differences are subtle, however, and not significant at upper/lower quartile levels. The null hypothesis 1 presented above cannot be rejected by these data due to the lack of systematic variation.
That LPI did not clearly separate 'not observed' from 'occurrence' of ejected sediment, or the 'minor' and 'moderate-severe' classes of severity, may simply highlight that such geotechnical parameters tend to be probabilistic (rather than deterministic) indices with a variety of possible shortcomings (e.g. Cubrinovski et al., 2019). Alternatively, other causal factors could have been involved in the severity of liquefaction manifestation and degree of damage. Trends can also sometimes become clearer by combining observations and calculations from multiple earthquake events (e.g. Maurer et al., 2014), rather than separating individual earthquakes as done here. Detailed examination of relationships between geotechnical parameters and liquefaction damage during different Canterbury earthquakes is provided elsewhere (e.g. Tonkin and Taylor Ltd, 2013; van Ballegooy et al., 2014bvan Ballegooy et al., , 2015avan Ballegooy et al., , 2015bCubrinovski et al., 2019).
We then analyzed the dynamic head change and total coseismic head of groundwater in deep aquifers against the response variables of observed ejected sediment within the categories of low (LPI ≤ 5), medium (5 < LPI ≤ 15), and high (LPI > 15) liquefaction vulnerability, using classes derived locally for Christchurch (Fig. 7). The approach uses LPI to isolate the combined variance in shaking, water table depth and soil profiles. Essentially it evaluates a null hypothesis 2 "that observed occurrence or severity of ejected sediment is unrelated to aquifer pressure when controlled for liquefaction vulnerability (based on combined shaking intensity, soil strength profile & saturation)". The dynamic head increases in the urban study area caused by the proximal Mw6.2 earthquake were mostly higher, and far more variable, than changes generated during the more-distal Mw7.1 earthquake (Figs. 1A, 7A cf. 7B). Although there are some spatial correlations at a regional scale (Fig. S2B, S2C), the variation in dynamic head change across the smaller urban area during the Mw7.1 was insufficient to separate populations of sites where ejected sediment was 'not observed' from those sites with 'minor' or 'moderatesevere' occurrences. The null hypothesis 2 cannot be rejected by data for the Mw7.1 earthquake. However, during the Mw6.2 earthquake, regardless of whether the liquefaction vulnerability was low, medium or high, the 'minor' and 'moderate-severe' liquefaction occurred where there was greater dynamic pressure in aquifers at depth. Median values for 'minor' and 'moderate-severe' response variables are nearly all above the 'not observed' upper quartiles, and 'moderate-severe' liquefaction mostly occurred where there was higher median head changes than where liquefaction was 'minor'. These dynamic head changes appear to carry through into the total coseismic hydraulic head condition, with differences at quartile-level significance in the populations of places where liquefaction occurred versus where it was 'not observed' during the Mw6.2 earthquake. Those sites where 'minor' or 'moderatesevere' amounts of sediment were ejected also had higher total coseismic head (in these instances where the total head exceeded 5 m). During either the Mw7.1 or Mw6.2 earthquake, moderate-severe sediment ejection was mostly restricted to places where the total coseismic pressure was above ground (10th percentiles >0, Figs. 7C, D). Although the null hypothesis 2 cannot be rejected by Mw7.1 data, the alternate "that observed occurrence or severity of ejected sediment is correlated to aquifer pressure" can be accepted for the Mw6.2 earthquake at quartile-level of significance, regardless of whether sites had low, medium or high liquefaction vulnerability. Trends in Fig. 7B & 7D suggest the LPI values may be missing an important factor that influences the occurrence of liquefaction. Given the assessment is non-spatial, perhaps the wider spatial correlation described in Section 3.1 (Fig. 4) is not merely coincidental.

The role of soil behavior
Differences between vulnerability predicted from geotechnical indices and observed manifestation of ejected sediment and/or liquefaction-associated land damage observed in Christchurch have been highlighted by several studies (e.g. Tonkin and Taylor Ltd., 2013;Maurer et al., 2014;van Ballegooy et al., 2014b). Geotechnical indices tended to overpredict damage in southwestern Christchurch (suburbs of Halswell, Hoon Hay, & Spreydon) and underpredict the severity observed in the northern township of Kaiapoi and some suburbs of eastern Christchurch, but over-and underpredictions are not always consistent from one earthquake to the next (Maurer et al., 2014;van Ballegooy et al., 2015b).
There are many potential reasons for this difference, that include: local site variability of geological and soil-strength parameters; anthropogenic modification or infrastructure; presence of silty layering and zones of partial saturation; uncertainty in water levels during the earthquakes; geotechnical testing after the earthquake, when water levels may have been different; and uncertainty in horizontal peak ground acceleration and calculations. It is also possible that vertical acceleration, surface waves or dissipated energy may cause fluid pressure increases that are not fully considered in vulnerability calculations (Davis and Berrill, 2001;Holzer and Youd, 2007). The standard triggering parameters are also simplified, such that interactions between different layers in the dynamic response, or pore-water pressure redistribution and water flow, are generally ignored (Cubrinovski et al., 2019). Differences can also be caused by the process and timing of liquefaction mapping. It may be incomplete, for example when some immediate clean-up has occurred before it was mapped, or it may be complicated by places where ejected sediment and water flow carried material laterally away from eruptive centers (Movie S1).
Some factors that could cause disagreement can be expected to vary from one earthquake to the next at each site, such as seismic loading and/or its characterization. Other parameters such as soil strength are more likely to be constant at a site, regardless of the earthquake. In a seminal study, Maurer et al. (2014) examined ~1200 sites and identified a trend of overpredictions for Christchurch soils with a behavior index Ic ≥ 2.2. They suggested that plastic soils with elevated silt-or clayfractions can limit flow of water to the surface because they are more resistant to piping and hydraulic fracture, resulting in local overprediction. Maps of geomorphological environment (Fig.8A, after Begg et al., 2015) and interpolated Ic10 (Fig.8B) both show demarcation between fluvial and lacustrine depositional environments in the urban study area: sandy-silt or clayey soil types are present in the west (Ic10 ≥ 2.2), whereas silty-sand or clean-sand behavior-types (Ic10 < 2.2) dominate coastal dunes and shorefront sites in the east.
Motivated by the need to understand the possible influence of soil behavior and aquifer pressure redistribution, we divided the geotechnical datasets into sites of underprediction where liquefaction vulnerability parameters were low (LPI ≤ 5) but ejected sediment was 'moderate-severe', and overprediction where vulnerability was calculated to be high (LPI > 15) but ejected sediment was 'not observed'. Of the 11,603 sites for which LPI is calculated for both the Mw7.1 and Mw6.2 earthquakes, around 14% had one incidence of either over-or under-prediction, but only 89 sites (0.8%) were consistently overpredicted and 36 (0.3%) underpredicted for both events. A map of sites with consistent underprediction (red circles) and overprediction (blue circles) shows relatively constant behavior over large areas (Fig.8B), whereas mapping of single-earthquake over-and under-prediction displays a similar pattern but with variability at the kilometer-scale (Fig. 8D). Populations of consistently underpredicted (n = 36) and overpredicted (n = 89) sites are clearly separated in terms of both Ic10 soil behavior and total coseismic head (Fig. 8C,F), with underprediction at places with lower Ic10 (inferred to be clean sands or silty-sands), or where there is higher aquifer pressure. Most of the consistently underpredicted sites are where deep aquifers have coseismic hydraulic head modelled to have been >6 m during the Mw6.2 earthquake (Fig.8E), and many of the smaller km-scale variations in over-and underprediction appear to coincide with changes in aquitard thickness (Fig.8E) and/or sedimentary facies (Fig.8A). But both trends have potential to contain artefacts of spatial association caused by southeast-northwest changes in coseismic aquifer pressure that are broadly coincident with east-west gradients in shallow soil-type (Fig. 8AB cf. Fig. 3EF) and/or seismic loading (Fig. S1).
To further examine the effects of soil behavior, we analyzed dynamic and total coseismic pressure in aquifers against observed ejected sediment ('not observed', 'minor', 'moderate-severe') during the Mw6.2 earthquake, with sites categorized by soil behavior index (Ic10 ≤ 1.8, 1.8 < Ic10 ≤ 2.2, Ic10 > 2.2) (Fig. 9A,B; see also Fig. S3C). Populations of site values where liquefaction was not observed have progressively lower median values of dynamic head change and total coseismic head relative to sites where liquefaction was not observed for each Ic10 bin. This may reflect the spatial coincidence, mentioned above, but a lack of differences in the medians for moderate-severe or minor liquefaction across the soil behavior categories is noteworthy. Perhaps more importantly, within each category of soil behavior there is a clear tendency towards higher median dynamic head change, or total coseismic hydraulic head, where ejected sediment was observed. The trends displayed in Fig. 9 suggest that although soil behavior is a locally important variable to incorporate in the prediction of liquefaction, larger-scale variations in aquifer pressure may also have been important for liquefaction during the Mw6.2 earthquake.
As a final test, we re-examined hypothesis 2 using a subset of sites with Ic10 ≤ 1.8, inferred to be clean sand sites, typically occurring towards the eastern side of the city (Fig. 8B) where soil profiles are unlikely to contain much interlayered silt, clay or peat. This controls any variance in soil behavior-type and reduces the potential for artefacts from coincidental spatial association. Dynamic head changes and total coseismic head in aquifers at these sites were assessed against ejected sediment during the Mw6.2 earthquake, categorized by liquefaction vulnerability (Fig. 10). Vulnerability was calculated to be low (LPI ≤ 5) for nearly one-quarter (1136/4944) of the Ic10 ≤ 1.8 sites during this earthquake. Where standard interpretation would suggest there should have been little or no liquefaction, 42% experienced 'minor' and 22% 'moderate-severe' liquefaction manifestation, and these were clearly where the aquifers below had higher dynamic pressure than where liquefaction was 'not observed'. The dynamic head changes (Fig.10A) correlate with the total coseismic hydraulic head, which show a similar pattern (Fig.10B). Null hypothesis 2 is clearly rejected by low (LPI ≤ 5) vulnerability data for the Ic10 ≤ 1.8 subset. As any coincidental spatial association should have been removed for this subset, the aquifer pressure does appear to directly affect both occurrence or severity of ejected sediment in this situation. The results also suggest that failing to take Fig. 9. Dynamic and total coseismic aquifer pressures at geotechnical testing sites within the Christchurch urban area during the Mw6.2 earthquake. Hydraulic head parameters are plotted against ejected sediment classes, factored against site soil behavior indices (average Ic10). Box and whiskers follow the same format and colour scheme as previous plots. account of confined aquifer pressure can result in local underprediction of liquefaction.

Background context
Earthquake-induced changes in permeability or rapid, essentially catastrophic, breaching of hydraulic barriers and seals is quite widely recognized in geological and geophysical literature (see Wang and Manga, 2010;Manga et al., 2012;Ingebritson and Manga, 2019;and refs therein). At deep crustal levels, perturbation of hydrologic systems is intimately linked and a controlling factor in earthquake processes, stress cycling, crustal strength and mineral deposit formation (e.g. Sibson, 1994;Ellsworth, 2013). Nearer the surface, earthquake perturbation of groundwater systems is regularly recorded in wells and monitoring networks (e.g. Elkhoury et al., 2006;Weaver et al., 2020) providing evidence of aquitard breach and/or switch of aquifers from confined to semi-confined Shi and Wang, 2016). Flow and changes in groundwater resource volumes depend on re-equilibration of local hydraulic gradients, which controls whether wells and aquifers will record long-term and sustained rises or falls (e.g. Rutter et al., 2016;Hosono et al., 2019).
But although earthquake hydrologic responses are widely observed, their effects are less-commonly of direct concern to engineering. Impacts are clearly important and can be hazardous in underground mines, nuclear waste-repositories, and groundwater supplies (e.g. Shimizu et al., 1996;Yoshida et al., 2000;Li et al., 2007). In the near-surface, redistribution of pore-fluid pressure can affect landslide stability (O'Brien et al., 2016) and was proposed as a mechanism contributing to triggering of liquefaction at far-field distances (Wang, 2007). There are numerous examples where earthquakes alter spring discharge and chemistry, or change stream flows, but these effects are commonly ephemeral when compared with the impact of rainfall-or groundwater-related flooding events (e.g. MacDonald et al., 2012). Substantive groundwater discharge at the Earth's surface requires positive vertical gradients and/or sustained hydraulic head at above ground pressure in permeable material, which can also be generated by dynamic or static stresses (Bonini, 2020). Eruptive 'mud volcano' slurries of fluidized sediment and gas (e.g. Bonini et al., 2016;Kassi et al., 2017) are clear examples where hydrologic responses to earthquakes (as well as anthropogenic-or spontaneous-triggering) can translate directly into local hazards and surface changes with impacts for engineering. Some of the more-vulnerable hydrogeological settings are arguably where permeable aquifers are interbedded with fine-grained confining and semi-confining aquitards along coastal margins and in rift valleys (e.g. Yechieli and Bein, 2002). Here aquifers are recharged at higher elevations and sit in close proximity to active faults, and unless there has been appreciable groundwater extraction, have potential to cause impact where hydraulic heads are near or above the land surface. With positive hydraulic gradients and head in excess of 5 m in aquifers beneath Christchurch (Fig.2E), the Canterbury plains and coastal aquifer system clearly meet these criteria.

Mechanistic interpretation
The shallow aquitards in eastern Christchurch comprise interbedded Upper Springston Formation (sandy gravel channel deposits and overbank silts) and Christchurch Formation (marine/estuarine peat, clay, silt and sand) (Brown and Weeber, 1992;Begg et al., 2015). The aquitards are < 40 m thick beneath the eastern suburbs (Talbot et al., 1986;Weeber, 2008), where they confine positive artesian pressures of at least 2 m, and locally > 5 m (Fig. 2). The footage in Movie S1 (Supplementary Materials) highlights some major eruptive centers in this area that were spaced around 50 m apart or more, reaching peak flows delayed by around 2-5 min after shaking started. In many instances flows continued for hours after the event and well beyond the length of posted videos. These examples demonstrate the large volumes (~10 5 L) of water Fig. 10. Dynamic and total coseismic aquifer pressures during the Mw6.2 earthquake at a subset of geotechnical testing sites where the Ic10 ≤ 1.8 (soil behavior-type approximating clean sand). Hydraulic head parameters are plotted against ejected sediment classes, factored against liquefaction vulnerability categories defined by LPI values. Box and whiskers follow the same format and colour scheme as previous plots.
issuing from eruptive centers, carrying sediment, with peak flows locally exceeding 10 Ls − 1 . We find it difficult to envisage how lateral flow and expulsion from within the shallowest confining aquitard units, themselves locally liquefied, could remain isolated from dynamically pressured artesian aquifers just some tens of metres below. Is it possible that a proportion of high-discharge groundwater flows and surface flooding could be the result of physical leakage from the deeper aquifer system, with non-Darcian flow carrying liquefied sediment and/or products of subterranean erosion to the surface slightly delayed after shaking?
Although there are clear differences in the amount of sediment manifested across the study area during the Mw6.2 Christchurch earthquake compared with the Mw7.1 Darfield earthquake, and differences in the seismic loading magnitude and gradients across the city (Fig.S1), the overall shape of spatial probability distributions derived for ejected sediment have many similarities between the earthquakes (Figs. 5 cf. S2). We confirmed expectations from engineering (Ishihara, 1985;Idriss and Boulanger, 2008) that sediment ejection and surface manifestation occurred predominantly where the shallow unconfined water table was near-surface, with very few occurrences where it was > 4 m deep, but found water table depth had little bearing on the severity sub-classes of minor or moderate-severe sediment ejected in the Christchurch urban area (Fig. 5A). Mapped sediment occurrence has a strong spatial correlation with dynamic head change induced by the Mw6.2 earthquake (Figs. 4A, 5B, 7B), which locally added up to 5.6 m to the total hydraulic head, but correlation was not as clear with < 2 m changes and smaller amounts of liquefaction induced by the Mw7.1 event (Fig. S2B, Fig. 7A). Although ejection only occurred where dynamic head changes were positive, the shaking-induced changes (increases) were relatively constant across the urban region during the Mw7.1 earthquake and do not appear to have a first-order control on any spatial variation in the sediment ejected during the first, more-distal, earthquake.
Relationships between ejected sediment and total coseismic head in deeper aquifers provide notable correlations in both occurrence and severity of ejected sediment, which rose sharply where total coseismic head became positive (artesian) relative to ground level (Figs. 5C, S2C). The spatial probability distributions for the Mw6.2 and Mw7.1 earthquakes were similar in their shape, but the variation across the study area caused by the more-distal Mw7.1 earthquake was much less pronounced. Manifestation during the Mw6.2 earthquake appears related to the dynamic head change and total coseismic head in groundwaterregardless of whether geotechnical parameters (accounting for effects of shaking, water table depth and soil profile strength) indicated sites to have low, medium or high liquefaction vulnerability (Fig. 7B, D). Although initial nonspatial treatment of the entire geotechnical dataset potentially still contains some covariate associations that are spatial artefacts, the examination of ordinal multinomial responses using selected groups and subsets (Fig.10) maximized the opportunity to elucidate direct correlations and causal relationships. The nonspatial treatment of data supports the alternate hypothesis 2 that observed occurrence or severity of ejected sediment was locally related to aquifer pressure in some circumstances. Sites that exhibit clean sand-type behavior (Ic10 ≤ 1.8) and should have low vulnerability (LPI ≤ 5) are commonly underpredicted. A mechanistic explanation is that positive hydraulic head relative to ground, not considered within the LPI calculation, provides a driving force enabling groundwater and entrained sediment to be carried to the surface.
Positive spatial correlations between confining sediment thickness and ejected sediment occurrence (Fig. 5E) may simply confirm an expectation that thick sediment profiles are prone to liquefaction and/or producing ejecta. However other studies have observed and modelled significant changes in hydraulic conductivity (Dudley Ward, 2015;Weaver et al., 2019), upward passage of groundwater (Cox et al., 2012;Gulley et al., 2013: Dudley Ward, 2015, and long-term offset in confined and unconfined groundwater levels (van Ballegooy et al., 2014a;Quigley et al., 2016;Rutter et al., 2016) induced by the Canterbury earthquakes. If the aquitard initially impeded groundwater leakage it may also have aided buildup of aquifer pressure before the onset of liquefaction. As cyclic loading continued it could have then increased its hydraulic conductivity as it fractured or layers within it liquefied, then contributed abundant silty-sand to be entrained and carried by flow to the surface. The extent and location of sediment could then directly correlate to levels of pressure in the aquifer relative to ground. This would explain both the correlation between confining sediment thickness and ejected sediment occurrence, and why sites where liquefaction was underpredicted by geotechnical parameters had median coseismic groundwater pressures that were high relative to sites where liquefaction was overpredicted (Figs. 8,9).
Based on correlations with hydrogeological parameters and observations during the Mw7.1 and Mw6.2 earthquakes (Figs. 5, S2, S3; Movie S1), we propose that ejection of groundwater and sediment in Christchurch, and quite a lot of material referred to as 'liquefaction', involved both pore water expelled by shallow volumetric densification and a component of deeper groundwater that was unable to remain confined in aquifers during these ground-damaging earthquakes. We envisage confinement that was compromised by subsurface fracturing and/or liquefaction (cf. Shi and Wang, 2016;Wang et al., 2016), such that groundwater escaped from aquifers dynamically pressured by the earthquakes. Underground erosion, suffusion and piping along fractures and liquefied horizons entrained sediment in groundwater en-route to the surface, creating subsurface voids that contributed to differential settlement. Although clearly not a prerequisite for liquefaction ejection, which has been observed worldwide in the absence of confined groundwater (e.g. National Research Council, 1985;Terzaghi et al., 1996;Youd et al., 2001), such a process explains our wider-scale spatial observations that surface manifestation and ejection of sediment and water: 1) was enhanced where total coseismic groundwater pressure in deep aquifers was above ground (artesian); 2) nearly all occurred where the coseismic vertical head gradients were upwards (positive); 3) had occurrence and consequence that are locally difficult to predict, potentially explaining kilometer-scale disparities in the absence of other local influences (e.g. such as silty-soil interlayering or partial saturation); and 4) was relatively extreme during the Mw6.2 by most worldwide comparisons.

Application of results
On the basis of well-known laboratory experiments, theory, and observations of ground surface subsidence measured by LiDAR surveys and tectonic slip models (e.g. Ishihara, 1985;National Research Council, 1985;Terzaghi et al., 1996;Youd et al., 2001;Idriss and Boulanger, 2008;van Ballegooy et al., 2014b;Quigley et al., 2016;), it is generally concluded that manifestation of sediment in Christchurch was simply the direct consequence of liquefaction, in which the cyclic strain induced by earthquake shaking in loose to medium dense sand resulted in an increase in pore water pressure and volumetric densification that squeezed out excess pore water and ejected water and sediment at the ground surface. Whether the process of liquefaction involved any distributed release of aquifer pressure and Darcian flow, as proposed by Wang (2007), remains an open question that cannot be answered with groundwater monitoring data presently available for Canterbury (see Supplementary Materials: Discussion). The term 'liquefaction' became a New Zealand household noun to describe any unwanted water and sediment that repeatedly inundated people's homes and suburbs. It is, however, quite difficult to establish how widespread volumetric densification was and the length-scales at which it occurred. Differential LiDAR surveys provide a measure of the land's total vertical elevation change (e.g. Hughes et al., 2015), but in Christchurch a significant and imprecise component of tectonic motion (uplift and subsidence on blind faults) needs to be distinguished from more-localized sediment compaction before the depths and length-scales of densification can estimated (e.g. Beavan et al., 2011;Tonkin and Taylor Ltd, 2013). There has also been little evidence found for changes of geotechnical ground strength properties which correlate with density at sites that have been repeatedly analyzed (e.g. Lees et al., 2015;Orense et al., 2012;Russell et al., 2015 Appendix E).
We found the spatial probabilities of ejected sediment increase with Liquefaction Potential Index (LPI) as expected (Figs. 6, S3). However, with widespread differences between predicted vulnerability and observed liquefaction damage (van Ballegooy et al., 2015b), the LPI spatial probability distributions did not seem quite as pronounced as those calculated for hydrological parameters (Fig. 5). Cumulative curves of ejected sediment 'occurrence' versus 'not observed' for the geotechnical indices seem less-separated than for hydrological parameters, due in part to the long tail nature of their probability distributions. Through nonspatial assessment using ordinal multinomial responses, it was not possible to reject a null hypothesis 1 "that occurrence or severity of ejected sediment is unrelated to liquefaction parameters when controlled for aquifer pressure" for either earthquake. Instead, differences in ejected sediment are associated with both dynamic pressure changes and total coseismic head during the Mw6.2 earthquake, regardless of liquefaction vulnerability (Fig. 7), which appear to be causative at low vulnerability (LPI ≤ 5) sites with sand-type behavior (Ic10 ≤ 1.8) (Fig. 10). An alternate hypothesis 2 "that observed occurrence or severity of ejected sediment is related to aquifer pressure" was supported by data for the proximal Mw6.2 event. So, while the current geotechnical indices appear to provide a method of evaluating site liquefaction probabilistically, and the derived parameters are widely-adopted for defining general areas with soil conditions of concern, the hydrological parameters for Christchurch can indicate where sediment manifestation may be exacerbated and/or underpredicted. The hydrological parameters can potentially be used for defining suburb-scale areas of elevated hazard.
Regardless of whether distributed (coseismic) or discrete (delayed) passage of groundwater occurred, the positive spatial correlation observed between ejected sediment occurrence and confining sediment thickness (Fig. 5E) will be important to hydrogeologists, as it suggests the shallow aquitard in Christchurch may no longer have provided protection from geopressures generated when the aquifer system was impacted by a large earthquake nearby. Although not all areas with artesian pressure leaked, and we are not yet able to predict where leakage might be expected at exact land-parcel scale, we suggest leakage areas likely experienced a cascade of non-linear erosive processes that temporarily transformed the shallow sediment structure. Defining areas of above-ground (artesian) groundwater pressure, as both static and potential coseismic dynamic head, provides a first-pass suburb-scale (> 1 km) indication as to where such processes could occur and exacerbate the liquefaction-induced ground damage. There is potential for mapping and pre-assessing whether geotechnical indices of liquefaction potential will over-or under-predict damage (cf. Fig. 8).
Unprecedented groundwater monitoring and geotechnical datasets in Canterbury have enabled demonstration that the extent and severity of sediment ejected to the surface in Christchurch spatially coincided with deep aquifer pressure and, in particular, the places where groundwater was confined at artesian (above ground) pressure. Although it is difficult to quantify the contributions of aquifer leakage cf. liquefaction per se to damage, they involve distinct mechanistic processes that operate together, are likely non-linear and sporadic, but have potential to be mitigated separately. Confined artesian groundwater presents a unique situation where sand, silt and dirty water can potentially escape to the surface following ground damaging levels of shaking, then drive repeated inundation and subsurface erosion in subsequent aftershocks once pressures have recovered. As such, where there is sufficiently strong-shaking of loose soils, artesian aquifer systems may host a distinct, and as yet little recognized, additional geological hazard. Christchurch was one of the most pervasive and extensive urban liquefaction events ever experienced, that has been captured by exceptional groundwater and geotechnical datasets. Extracting information on processes, however, is not entirely straightforward, because both known and proposed contributory factors covary in space. A logical step to extend this research will be to look more widely and critically at other examples worldwide, for instances where liquefaction occurrence or damage associated with ejection may also have been exacerbated by a hydrogeological pre-condition, but where spatial associations and hydrogeology may be different. Laboratories might also investigate effects when a localized component of basal fluid pressure is added to liquefaction experiments.
A corollary is that during the evaluation of liquefaction potential, the wider hydrogeological environment does need to be considered. Built over aquifers with substantial artesian groundwater pressures that were dynamically increased by the earthquakes, the severity of the observed liquefaction-induced ground damage in Christchurch during the 22 Feb 2011 Mw6.2 seems to have been an extreme example that may not be widely applicable elsewhere. Not all confined systems will have the same sensitivity given similar levels of ground shaking. Sustainable development on artesian aquifer systems will be dependent on: the thickness of the non-liquefying surface layers (often influenced by the depth to the shallow groundwater surface); the grainsize, sorting and strength characteristics of liquefying soil layers; the thickness and integrity of the confining layers hydraulically isolating the artesian aquifer; and the pressure and availablestorage characteristics of the groundwater confined in the aquifer below. There should also be expectation that hazards may change seasonally or with sea-level rise, however in many urban environments they may have already been unconsciously mitigated by groundwater abstraction and aquifer drawdown.

Conclusions
A variety of hydrological, ground damage and geotechnical data were compiled for Canterbury spanning an area that includes both inland unconfined/semi-confined and coastal confined aquifers. A network of groundwater monitoring enabled definition of preearthquake static groundwater levels, dynamic changes induced by earthquakes, and total coseismic hydraulic head in the aquifers during two major Mw7.1 and Mw6.2 events of the 2010-2011 Canterbury earthquake sequence. In a wider spatial analysis of data, we presented evidence that the spatial probability of sediment ejection associated with liquefaction during the earthquakes was greatest where absolute groundwater pressures in the deep aquifers reached levels that were above ground and vertical head gradients were positive.
As spatial correlations have potential to be coincidental associations rather than reflecting direct causal relationships, we then undertook a nonspatial assessment of ordinal multinomial responses in data from sites where geotechnical testing has been completed in the Christchurch urban area. Although spanning a smaller range of hydrogeological conditions and biased towards places that experienced liquefaction damage, by grouping geotechnical data based on liquefaction vulnerability indices, we isolated the combined effects of variations in shaking, water table depth and/or soil profile strength. For the Mw7.1 Darfield earthquake results were equivocal, possibly because the dynamic pressure changes in aquifers were smaller and showed little spatial variation across the restricted urban study area. It was clear, however, that during the Mw6.2 Christchurch earthquake there was significantly higher coseismic hydraulic head in aquifers beneath places where ejected sediment was 'moderate-severe' or 'minor', compared with 'not observed'. Regardless of whether the ground had low-, medium-or high liquefaction vulnerability, there are differences in groundwater pressure associated with occurrence and severity of ejected sediment are consistent during the Mw6.2. A subset of sites with clean sand-type behavior (Ic10 ≤ 1.8) and low predicted vulnerability (LPI ≤ 5) show relationships that are unequivocal and very unlikely to be a coincidental spatial association. In specific situations it seems clear that aquifer pressure can directly affect occurrence of ejected sediment and contribute to local underprediction of liquefaction vulnerability.
We collated evidence and argue there are logical mechanistic reasons for the observed relationships and that potentiometric levels provided the driving force for open-system loss of groundwater, leading to subterranean erosion, suffusion and piping. Leakage from artesian aquifers is certainly not a prerequisite for the surface ejection of sediment and water, as pore-pressures can be increased by cyclic shearing of loose to medium dense sandy soils in the total absence of confined groundwater. But where shaking is sufficiently strong for aquitards to be compromised and vertical flows to initiate, we envisage a strongly non-linear cascade of subsurface erosion processes through a soft sediment pile. The release of artesian water may in turn exacerbate any near-surface liquefactionrelated damage. The 22 Feb 2011 Christchurch Mw6.2 earthquake is an example that may well have been extreme but provides rationale for reevaluating other examples of liquefaction worldwide and consideration of ground performance in relation to wider hydrogeological conditions and setting.

Data availability
An ArcGIS 10.5.1 geodatabase and Microsoft Excel files containing data related to this article can be found at, Zenodo.com .

Statement of contributions
The authors do not have any competing financial interests. The study was designed by S.C.C., S.vB., M.M. and H.J.R. following an initial hypothesis by our colleague Chi-Yuen Wang (University of California Berkeley). Water level, shaking and liquefaction triggering parameters were calculated by H.J.R., C.H. and V.L., respectively. Modelling of groundwater pressures, spatial probability distributions and statistical analysis was completed by S.C.C., D.S.H. and A.K.G., with various technical support. The principal author of the manuscript was S.C.C., but all coauthors discussed the content, provided comment, and support the results. Correspondence and requests for materials can be sent to S.C.C.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.