Seismological and Hydrogeological Controls on New Zealand-Wide Groundwater Level Changes Induced by the 2016 Mw 7 . 8 Kaik ō ura Earthquake

The 2016 Mw 7.8 Kaikōura earthquake induced groundwater level changes throughout New Zealand. Water level changes were recorded at 433 sites in compositionally diverse, young, shallow aquifers, at distances of between 4 and 850 km from the earthquake epicentre. Water level changes are inconsistent with static stress changes but do correlate with peak ground acceleration (PGA). At PGAs exceeding ~2m/s, water level changes were predominantly persistent increases. At lower PGAs, there were approximately equal numbers of persistent water level increases and decreases. Shear-induced consolidation is interpreted to be the predominant mechanism causing groundwater changes at accelerations exceeding ~2m/s, whereas permeability enhancement is interpreted to predominate at lower levels of ground acceleration. Water level changes occur more frequently north of the epicentre, as a result of the fault’s northward rupture and resulting directivity effects. Local hydrogeological conditions also contributed to the observed responses, with larger water level changes occurring in deeper wells and in well-consolidated rocks at equivalent PGA levels.


Introduction
Central New Zealand has experienced several moment magnitude (M w ) 7 or larger earthquake ruptures in the last 200 years, the largest of which were the 1848 M w 7.4−7.7 Awatere earthquake [1,2], the 1855 M w 8.1+ Wairarapa earthquake [3], and the 1888 M w 7−7.3 Hope earthquake [4].The most recent M w 7+ event in the region was the M w 7.8 Kaikōura earthquake, which occurred on 14 November 2016 at 00:02:56 (New Zealand Daylight Time, NZDT).The earthquake was the largest in almost a decade of seismic disruption in New Zealand [5], following four decades of relative quiescence [6,7].The Kaikōura earthquake followed the 2010-2011 Canterbury earthquake sequence, which included the devastating M w 6.2 aftershock of 22 February 2011 [8,9], and the Cook Strait earthquakes of 2013 [10,11].
The Kaikōura earthquake initiated in the Waiau Plains in northern Canterbury, New Zealand, on an oblique thrust at a depth of ~15 km.It propagated upwards and northwards, forming a complex rupture along a ~180 km zone that involved at least 21 faults [12,13], producing widespread shaking reaching maximum Modified Mercalli Intensity of IX.The maximum horizontal surface displacement was 10 m in a dextral sense [14,15] and vertical displacements were highly variable, causing widespread coastal uplift/subsidence between +6.5 and −2.5 m [16].The rupture took over two minutes [5], and ground acceleration exceeded ~1 g in the near-source region [17].
Earthquakes are known to induce groundwater level changes with varying amplitude, duration, and polarity of response (e.g., [18][19][20]).Earthquake-induced static and dynamic stresses (e.g., [21]) and local hydrogeological factors [22] influence groundwater level responses.In this first report on the Kaikōura earthquake's hydrogeological effects, we document an extensive dataset of groundwater level responses in order to assess the level to which earthquakedriven and local hydrogeological factors contributed to the national scale groundwater level response.This earthquake hydrology dataset is the most extensive compiled in New Zealand, surpassing that compiled by Cox et al. [23] following the M w 7.1 Darfield (Canterbury) earthquake, and is of comparable size to datasets compiled following the 1999 M w 7.5 Chi-Chi earthquake in Taiwan [24,25].In comparison to other case studies, the Kaikōura earthquake dataset is also significant as it includes groundwater level responses collated from a variety of hydrogeological settings.

New Zealand Hydrogeology
There are a wide variety of compositionally diverse aquifers in New Zealand that were formed in a range of depositional environments, mainly during the Tertiary and Quaternary geological periods.Sedimentary aquifers were formed from terrestrial (e.g., glacial, aeolian), shallow (e.g., deltaic, carbonate), and deep marine sedimentary deposits, with mixed clast composition reflecting the varied source rock geology [26].Transgression/regression sequences created overlaps of terrestrial and marine deposits in coastal areas [27], resulting in large variations in vertical and horizontal permeability and storage capacity [28].There are also volcanic aquifers of ignimbrite and scoriaceous deposits formed in the central North Island from Taupō Volcanic Zone caldera eruptions and basalt flows, respectively.Groundwater is also stored in Mesozoic-Cretaceous metamorphic bedrock that forms mountain ranges in both North and South Islands but is seldom utilised for water supplies due to the low permeability and yield of the bedrock [26].

Hydrological and Seismic Data
The hydrological monitoring infrastructure in New Zealand is extensive due to the country's heavy reliance on groundwater for agriculture [26].Regional councils manage and maintain networks of boreholes instrumented by data loggers and pressure sensors that measure groundwater levels, representing piezometric levels.The monitoring infrastructure is not evenly distributed as groundwater tends to be extracted from the most permeable areas, predominantly gravel aquifers, and pumping test data are not readily available.Where repeated pumping test data are available, the inferred transmissivities are ambiguous [29], highly variable [30], and catchment dependent [31].
We collected groundwater level time series data spanning the Kaikōura earthquake from 433 monitoring wells distributed throughout New Zealand (Figure 1).Groundwater levels are sampled at intervals ranging from one minute to three hours, with 15 minutes being the most common.Wells were classified on the basis of depth and local aquifer pressures, using the terms "water-table," "nonflowing," and "flowing artesian" based on the groundwater level relative to ground prior to the earthquake at each site.The majority of wells are less than 100 m deep, which are shallow in comparison with those considered in international studies (e.g., [32,33]).Screens typically range from 1 to 10 m in length, and well diameters generally vary from 50 to 300 mm, but overall, the construction method and well completion do not vary appreciably across the regions, resulting in similar style of earthquake-induced responses.The aquifers mostly occur in young sedimentary deposits that have varying levels of consolidation, with few aquifers in igneous or metamorphic rocks.A notable exception are the geoengineered schist landslides in Cromwell Gorge [34], which have fracture permeability and water tables that have been displaced from their natural condition by drainage tunnels.In contrast to these examples from New Zealand, many international studies of earthquake hydrology provide datasets on aquifers in well-consolidated and crystalline rocks (e.g., [35][36][37]).
We gathered rainfall and atmospheric pressure data from 65 climate stations, which are managed and maintained by the National Institute for Water and Atmospheric Research.Barometric corrections were applied where necessary using data from the nearest climate station, leaving monitoring well water level data relative to water level only.Rainfall data were qualitatively compared to groundwater level data before, during, and after the Kaikōura earthquake, to assess any effect on observed water level changes.
Earthquake-induced water level changes can be split into coseismic and postseismic components relative to the period of shaking at the monitoring site.The primary responses observed in this study were postseismic, as most of the records were collected at 15-minute intervals, the first occurring after the shaking had stopped.Water level changes were classified as either having no response, responding transiently and returning to preearthquake levels within two hours, or responding persistently (increasing or decreasing) and reequilibrating at a new postearthquake level lasting several days or longer (Figure 2).For each water level change, the polarity, amplitude, and duration until reequilibrated or returned to preearthquake levels were recorded.The amplitude of persistent water level changes was estimated from the difference between the reequilibrated and preearthquake groundwater level.At this longer time scale (than coseismic responses) and since most wells are specifically for monitoring (rather than production), we assume that any wellbore storage effect has already passed in water table wells, but some observed responses in confined aquifers could result from changes in transmissivity.Water level time series were removed from the analysis in cases in which earthquakeinduced damage to the well compromised or prevented the recording of water level.
The following general terms have been adopted for monitoring wells with respect to earthquakes [38]: "near field" denotes wells within one rupture fault length and "intermediate field" represents wells from one to ten rupture fault lengths away.Considering the Kaikōura earthquake rupture zone extends ~180 km northeast of the epicentre (Figure 1), 2 Geofluids the near field region extends ~360 km to the north and ~180 km to the south of the epicentre (Figure 2).We compiled all available New Zealand seismic recordings of the Kaikōura earthquake from 306 strong motion and broadband seismic stations.The seismic stations are part of the National Seismograph Network and Strong Motion Network (http://geonet.org.nz).Instrument responses were corrected, and a band-pass filter with transition bands of 0.10-0.25 Hz and 24.50-25.50Hz was applied.Shaking parameters, peak ground velocity (PGV) and acceleration (PGA), were calculated at each seismic station and interpolated to the monitoring wells using the nearest neighbour interpolation method [39].In the absence of borehole seismometers, we make an implicit assumption that surface shaking observed at the seismograph sites, interpolated locally as parameters at the well site, is applicable over the depth of the wells since they are mostly shallow.

Groundwater Level Changes Induced by the Kaikōura Earthquake
The Kaikōura earthquake induced groundwater level changes across New Zealand in the near and intermediate fields.
Water level changes were recorded at 433 sites from 4 to 850 km from the earthquake epicentre (Figure 2), at seismic energy densities [40] ranging between 10 -2 and 10 5 J/m 3 (Figure 3).Although the relationship of seismic energy density, magnitude, and epicentral distance was derived from earthquakes in southern California [41], the term is accepted internationally and used for seismo-hydrological phenomena (e.g., Weingarten and Ge, 2013; [36,42,43] 4, [44,45]).Earthquake-induced static stress changes imposed on the surrounding crust cause volumetric strain changes within aquifer systems.Studies have suggested that volumetric strain changes cause water level changes (e.g., [46][47][48][49]).The contractional and dilatational volumetric strain changes are inferred to increase and decrease water levels, respectively, with predictable amplitudes (e.g., [50]).However, several fault lengths away from the epicentre, water level changes are commonly larger in amplitude and inconsistent in water level change polarity (increase/decrease) with model predictions [38,51,52].
To calculate the mean static stress changes (σ kk ), we used the slip distribution of Clark et al. [16] to calculate the internal strain field based on the formulation of Okada [53] at 500 m depth across the epicentral region.Using a Poisson's ratio of 0.25 and a shear modulus of 30 GPa, we then convert the strain to stress at 500 m depth.Modelled fault slip geometries for the 21 ruptured faults were based on the surface rupture field [15] and geodetic data.Positive σ kk indicates tensional stress changes, and negative σ kk indicates compressional stress changes.
As expected, σ kk was most significant in the near field (Figure 4).In total, 386 wells were in areas of contraction, where σ kk ranged from −0.2 to −18,800 kPa.47 wells were in areas of dilatation, where σ kk ranged from 0.7 to 120 kPa (Figure 5).Static stress-induced compressional changes larger than 10 1 kPa occurred predominantly with water level increases (n = 98), and static stress-induced tensional changes between 10 1 and 10 2 kPa occurred generally with water level decreases (n = 13).Compressional changes smaller than 10 1 kPa occurred more frequently with decreases in water level (n = 53) than increases (n = 39), although the larger amplitudes were mostly increases.At lower than ±10 2 kPa regardless of the sign of σ kk , some water levels did not respond.Significant (>1 m) water level changes were observed in Cromwell Gorge, where σ kk was less than −10 0 kPa (Figure 5).5.2.Earthquake-Induced Dynamic Shaking.Dynamic stress changes caused by the passage of seismic waves are shortlived and decay at ~1/r 5/3 , meaning that they can be significant from the near to far field [44,45].Dynamic stressinduced deformation depends on loading rates and cycles of inertial forces [45] and is influenced by rupture directivity and radiation patterns [40].
5 Geofluids peak dynamic stress, which has been utilised in permeability enhancement experiments (e.g., Roberts, 2005; [56]).Dynamic stresses are thought to control the incidence of shear-induced consolidation and dilatation [57][58][59], vertical (e.g., [40,60,61]) and horizontal (e.g., [62]) enhancement of permeability.Peak dynamic stress (σ D , GPa) can be expressed [63] as in terms of the shear modulus (μ S , GPa), shear wave velocity at the monitoring well (v s , m/s), and maximum peak ground velocity (PGV, m/s).The maximum PGV of horizontal and vertical motions has been used here as a proxy for dynamic stress (equation ( 1), [62,64,65]).The assessment of liquefaction, an earthquake-induced hydrogeological response (e.g., [40]), is typically based on considerations of changes in stress and uses the horizontal (geometric mean) peak ground acceleration (PGA) [66,67].PGA also correlates with the amplitude and occurrence of groundwater level changes [68,69].We compared the horizontal (geometric mean) PGA to groundwater level changes.Comparing both PGV and PGA to groundwater level changes allowed an assessment of whether velocity or acceleration was the better indicator for determining the polarity, amplitude, and/or duration of water level change.
PGV and PGA were the highest in the near field, adjacent to the fault rupture, and notably high in areas of Marlborough, Wellington, Manawatu-Wanganui, Taranaki, and Hawkes Bay regions.The extent of the fault rupture to the north of the epicentre and the northward rupture directivity of the earthquake resulted in a slower decay of PGV and PGA north of the epicentre and a more rapid decay of PGV and PGA south of the epicentre.The higher PGV and PGA north of the epicentre correlate with more frequent larger 6 Geofluids amplitude water level changes.North of the epicentre, 60% of water levels changed persistently, while south of the epicentre 48% of water levels changed persistently (Figures 6-8).At monitoring wells, the PGV ranged from 0.002 to 0.9 m/s and the PGA ranged from 0.01 to 11.8 m/s 2 .Above a PGV of ~0.3 m/s and a PGA of ~2 m/s 2 , the majority of persistent water level changes were increases, and the median response reequilibration times were 1455 and 585 minutes, respectively.Below a PGV of ~0.3 m/s and a PGA of ~2 m/s 2 , there was no preferred polarity, and the median response reequilibration time was 75 minutes below both thresholds.PGA differentiates polarity behaviour more clearly than PGV, with 32 persistent water level increases above ~2 m/s 2 , compared to 11 persistent water level increases above ~0.3m/s.Some monitoring wells that experienced low levels of shaking showed a larger water level change than wells that experienced high levels of shaking, such as wells in Cromwell Gorge (Figure 7).This may be due to local hydrogeological conditions that partly contribute to the characteristics of water level changes.

The Influence of Local Hydrogeological Factors on Water Level Changes
Although distance from the epicentre provides a rough approximation of how a monitoring well may respond to an earthquake, some monitoring wells respond with larger amplitudes of water level change than others at the same epicentral distance.Aquifer properties are likely to partly control earthquake-induced water level changes.The higher the transmissivity of the formation, the more the well water level change will be reflective of the formation pressure change, with the time taken for water levels to reequilibrate or return to preearthquake levels being governed by flow properties [70].
Since PGA differentiates polarity behaviour more clearly than PGV (Figure 7) and that earthquake-induced dynamic shaking occurred at all monitoring wells, we scaled the absolute water level change amplitudes (WLC, m) to the horizontal (geometric mean) PGA (m/s 2 ) experienced at each monitoring well (WLC/PGA, s 2 ).To further understand the local hydrogeological factors in contributing to water level changes, we compared WLC/PGA to the depth of monitoring well and to the average shear wave velocity as a representation of the degree of confinement and expected dynamic rock strength behaviour at each site.
6.1.Depth.Earthquake-induced water level changes vary with the degree of confinement in aquifer formations.Water level changes in unconfined aquifers are generally smaller than those in confined/semiconfined aquifers [50], as a result of the specific yield of unconfined aquifers being higher than the storativity of confined aquifers [71].Changes in monitoring well water levels in unconfined aquifers reflect the (de)watering process of pore spaces above/below the water table, which forms the aquifer upper boundary.Changes are also influenced by well storage effects, the impact of which increases with decreasing transmissivity [50].Assuming undrained conditions [64], changes in confined aquifer pore pressure (p, GPa), given by p = −Bσ kk /3, are proportional to mean stress (σ kk , GPa), and the magnitude of pore pressure change is governed by Skempton's coefficient (B) [50].We used the depth of each monitoring well (4 m to 1.4 km) as a first-order proxy for the degree of confinement of the aquifers studied.The monitoring wells are generally shallow with the median depth being 24 m.At any depth, WLC/PGA typically varies by two orders of magnitude.Generally, the deeper the monitoring well, the more sensitive the  9 Geofluids well was to water level change.A similar observation was recorded by Rutter et al. [29], who found that water level changes occurred more consistently in wells of greater than 80 m depth, compared to shallower wells.At depths less than 10 metres, WLC/PGA broadly ranged from ~10 -3 to 10 -1 s 2 .At depths between 10 and 100 metres, WLC/PGA mostly varied from ~10 -2 to 10 0 s 2 .Deeper than 100 metres, WLC/PGA largely fluctuated from ~10 -1 to 10 1 s 2 .All types of water level changes were observed over the entire depth range as shown by the kernel density estimate.Monitoring wells in Cromwell Gorge geoengineered landslides were generally the most sensitive (Figure 9).The WLC/PGA for unconfined water table aquifers does not exceed 10 0 .6.2.Shear Wave Velocity.Seismic waves attenuate differently in different rock types [38], and it is therefore unsurprising that hydrogeological factors partly control earthquakeinduced groundwater level changes.Site effects and varying crustal structures have differing efficiencies of dynamic shaking [40].Saturated unconsolidated media damp highfrequency motions, while amplifying lower frequencies [25,72,73].In contrast, consolidated/crystalline aquifer material amplify high-frequency shaking more than low-frequency shaking [74].Sensitive hydrological sites screened in wellconsolidated or crystalline rocks can respond to teleseismic earthquakes over 1000 km away and typically do so with hydroseismograms (e.g., [32,36,62,75,76]).Where static and dynamic stresses are insignificant, hydrogeological factors must partly control the hydrological response.The geoengineered schist landslides of Cromwell Gorge respond hydrologically with large amplitude, although static and dynamic stresses are insignificant (Figure 8).Where static and dynamic stresses are significant, hydrogeological factors may influence hydrological responses but are probably less substantial.In the near field, responses are not entirely consistent with the magnitude of static or dynamic stresses (Figure 8).
To reflect the varied hydrogeological conditions of the aquifers studied, we used site average shear wave velocity (Vs 30 , m/s) over depths between 0 and 30 m [77] at each monitoring well.Shear wave velocities range from 40 to 1040 m/s, representing very soft soil to weak rock, respectively [78].Vs 30 is typically towards the lower end of this range, the median Vs 10 Geofluids typically ranged from ~10 -3 to 10 0 s 2 .When Vs 30 was close to 1000 m/s, WLC/PGA was mostly between ~10 -2 and 10 1 s 2 .
There is no correlation between polarity of water level change and Vs 30 at any depth.Monitoring wells in Cromwell Gorge were particularly sensitive in deep soil and weak rock site classifications.There are also numerous sites that did not respond to the earthquake at a variety of Vs 30 values (Figure 10).

Discussion
7.1.Earthquake-Induced Static Stress Changes.Static stress changes induced by the Kaikōura earthquake were only significant in the near field (Figure 8) and were variable in polarity and amplitude as a result of the complex rupture zone (Figure 4).Mean static stress change (σ kk ), used in this study, can be converted into volumetric strain change (ε kk ), ε kk = σ kk /E, with Young's modulus (E, GPa) of the lithology.For Young's modulus, an arithmetic mean of 13 GPa for sands and gravels was used [79].ε kk has the relationship with pore-pressure change (p) in a monitoring well, p = − 2GB/3 1 + v u / 1 − 2v u ε kk (e.g., [50]).To obtain a first-order estimate of the static strain changes and corresponding water level changes in this investigation, arithmetic means of shear modulus (G, 5.6 GPa), Skempton coefficient (B, 0.62), and undrained Poisson ratio (v u , 0.33) of sands and gravels were used [79], resulting in a coefficient 2GB/3 1 + v u / 1 − 2 v u value of ~9 GPa.Therefore, σ kk of ±10 0 and ±10 1 kPa predicts a ~10 cm and ~1 m water level change, respectively.However, most wells with σ kk of <±10 0 kPa had a larger water level change than predicted, even in cases of wells that had a consistent polarity with static stress changes.Only wells with σ kk larger than ±10 0 kPa had smaller water level changes than predicted, although the polarities were generally consistent with static stress changes exceeding ±10 1 kPa (Figure 5).
Considering water level change polarity or amplitude was inconsistent with static stress change predictions, it is unlikely that static stress changes contributed to water level changes in the near field.In the intermediate field, water level changes are larger than static stress change prediction, notably monitoring wells in Cromwell Gorge (Figure 5), which is consistent with other studies [38,45,52,80,81].These results are dissimilar to findings in other studies where static stress changes correlate with water level change amplitude and polarity [46][47][48][49][50]82]. High model uncertainty resulting from variable slip geometries and magnitudes between different models may partly account for poor fit between static stress and groundwater level changes, but we consider that groundwater level changes are unlikely to be caused by static stress changes.11 Geofluids greatly exceed a threshold of ~10 -4 [58], shear-induced dilatation occurs [83].Groundwater levels can decrease persistently as a result of the decrease in pore pressure and increase in porosity [84].Shear-induced dilatation may have occurred in close proximity to the Kaikōura earthquake fault rupture.However, as a result of the numerous landslides [85] and limited number of monitoring wells adjacent to the fault rupture, it is unclear to what extent shear-induced dilatation may or may not have occurred.
Shear-induced consolidation and liquefaction [40,83] occur at lower cyclic shear strains, still exceeding ~10 -4 [57,59].Shear-induced consolidation predicts that water levels increase persistently as a result of dynamic shaking, due to consolidating the sediment within an aquifer and decreasing porosity [83].The threshold at which persistent water level increases predominantly occurred was at a horizontal (geometric mean) peak ground acceleration (PGA) of ~2 m/s 2 and a maximum peak ground velocity (PGV) of 0.3 PGA differentiated persistent water level increases from other response types more clearly than PGV (Figure 7).Dynamic shaking also caused liquefaction to occur in the Marlborough [86,87] and Wellington [17,88,89] regions.These areas experienced a PGA exceeding ~2 m/s 2 and had predominantly persistent water level increases.
At low thresholds of dynamic shaking, earthquakeinduced flow velocities may be strong enough to dislodge colloidal particles [90].It has been postulated that dislodging of colloidal particles from preferential flow pathways may enhance horizontal permeability ( [62,84,91] [29]).The change in water level may originate in close proximity to a local pressure source, which could be produced by liquefaction [92], or elevated hydraulic heads [93].Earthquakeinduced flow velocities may be strong enough to dislodge colloidal particles [90].Laboratory experiments on pore unclogging [56,94] and groundwater colour changes [95] support permeability enhancement by dislodging colloids.The location of permeability enhancement could occur either up-or down-hydraulic head gradient of a monitoring well.If enhancement occurs up-hydraulic head gradient of a monitoring well, water level in the well would increase, as more flow is directed towards the well.If enhancement occurs down-hydraulic head gradient, water level in the well would decrease, as flow is directed away from the well.Therefore, if a sufficiently large number of observations are made, the permeability enhancement model would predict a statistically random occurrence in the polarity of the water level change [84].In this study, there was roughly an equal number of persistent water level increases and decreases at accelerations lower than PGA of ~2 m/s 2 , with water level change amplitudes generally less than 1 m.Furthermore, flowing artesian wells have positive responses in 33 out of 39 instances (see Supplementary Materials).This is consistent  [78].Site classifications are converted from Vs 30 [77].Monitoring well depths shallower than 30 m have larger opaque symbols as Vs 30 is representative.Monitoring well depths deeper than 30 m have smaller transparent symbols as Vs 30 is less representative.
with the model of enhanced permeability.Persistent water level changes may have occurred above a PGA of ~2 m/s 2 as a result of enhanced permeability but may be concealed by larger amplitude changes caused by shear-induced consolidation [38].
Water level changes caused by dislodging colloidal particles may be described as gradual sustained changes because the response reequilibration time to new postearthquake levels depends on the distance between the source of the blockage and the monitoring well [62,92].Water level changes caused by consolidation are more likely to be described as abrupt changes because immediate changes in pore pressure occur as a result of an undrained volumetric change around the well [38].In this study, the median persistent response reequilibration time for monitoring wells that experienced a PGA above ~2 m/s 2 (n = 33) was 585 minutes, while below a PGA of ~2 m/s 2 (n = 204) was 75 minutes.The two populations of response times were statistically different.Shear-induced consolidation may occur on a large scale (km) [83], as the threshold of shaking is likely to be exceeded over a wide area near the earthquake epicentre.However, permeability enhancement caused by dislodging of colloids may occur more frequently on a small scale (m) [62] in preferential flow pathways.The response time may reflect the scale at which these processes occur at but should not be considered definitive.
International case studies have shown that undrained consolidation/dilatation is the dominant mechanism in the near field, with enhanced permeability being dominant in the intermediate field [38].In the field, seismic energy densities larger than ~10 1 J/m 3 are above the threshold for undrained consolidation, with sensitive sites experiencing consolidation above ~10 -1 J/m 3 [38,96].For the Chi-Chi and Hengchun earthquakes, the transition from consolidation to enhanced permeability was inferred at ~10 1 J/m 3 [84].In this investigation, that equates to an epicentral distance of 100 km (near field), where water level changes were still predominantly increasing (Figure 8).The terms seismic energy density [40], near and intermediate fields [38], and one fault rupture length [97] do not take into account directivity of an earthquake.The occurrence of water level changes broadly correlated with the spatially asymmetric distribution of PGA (Figure 7), with persistent water level changes occurring 60% of the time north of the epicentre and 48% of the time south of the epicentre.Therefore, the directivity of an earthquake should be considered in the threshold for systematic comparison with other studies.As persistent water level increases predominantly occurred above ~2 m/s 2 (Figure 7), where liquefaction occurred, and the random polarity of responses below ~2 m/s 2 , we find that the transition between shear-induced consolidation and enhanced permeability may occur at a PGA of ~2 m/s 2 .This is not a definitive threshold and needs to be further compared to other earthquakes and monitoring sites before being adopted more widely.
At low levels of dynamic shaking, dynamic volumetric stresses cause pore spaces to dilate and compress, which can lead to transient pulses of pore pressure [92] and poroelastic deformation (e.g., [64]).In this case study, transient water level changes occurred mainly below a PGA of ~2 m/s 2 with small amplitudes of water level changes, (Figure 7), and are likely to be a result of poroelastic deformation (e.g., [64,92]).

Local Hydrogeological Factors.
Although shear-induced consolidation and enhancement of permeability are likely to have partly contributed to the polarity and occurrence of water level changes, the amplitude of water level change had a weak correlation with PGA (Figure 7).Some monitoring wells were more sensitive to water level changes than others at similar ground accelerations, which may be partially a result of hydrogeological factors.400 km south of the epicentre (Figure 8), monitoring wells in Cromwell Gorge screened in schist rock exhibited water level changes of ~±3 m.Yet, monitoring wells within 200 km of Cromwell Gorge, mainly screened in unconsolidated deposits, responded persistently with small amplitudes of less than 1 m.Cromwell Gorge monitoring wells are known to be sensitive to hydrological change from multiple earthquakes.Cromwell Gorge groundwater levels are depressed below equilibrium levels by pumping and gravity drainage due to infrastructure, and the sensitivity may in part reflect anthropogenic modification of the groundwater regime [34].Devils Hole, Nevada, is another example of a site that is sensitive to change, with responses occurring at low seismic energy densities of 10 -6 J/m 3 [75].Since there is a disparity between the global dataset [84] and case studies (Figure 3), it is possible that response thresholds vary widely as a result of hydrogeological factors.Different catchment responses may also affect the hydrological response to the earthquakes (e.g., [31]).
Depth is positively correlated with WLC/PGA (Figure 9), which concurs with studies that observed pronounced earthquake-induced water level changes in deeper confined aquifers compared to shallow unconfined aquifers [29,50,98,99].This is because the specific yield of unconfined aquifers is higher than the storativity of confined aquifers [71].Also, with depth, the transmissivity generally becomes lower and the strain sensitivity of the water head becomes larger [50].
WLC/PGA is positively correlated with Vs 30 (Figure 10), which is consistent with sites screened in well-consolidated or crystalline rocks responding hydrologically to teleseismic earthquakes [33,62,75,76].Therefore, regardless of the magnitude of the earthquake, it appears that monitoring wells may have a predisposition to show water level changes of certain amplitude relative to dynamic shaking, based on the strength of the rock they are screened in.This concept was adopted by O'Brien et al. [34] who concluded that aquifer systems have the ability to resist and recover from dynamic shaking which is consistent between earthquakes.WLC/PGA varies over several orders of magnitude at any given depth or Vs 30 .This large variation of WLC/PGA may result from variable permeabilities or well-aquifer coupling which may partly contribute to amplitudes of water level changes [34,45,90].
The fact that 127 monitoring wells that experienced different levels of shaking did not exhibit water level changes further demonstrates that hydrogeological factors do 13 Geofluids contribute to a monitoring well's capacity to exhibit a water level change.The well and monitoring system and/or the aquifer properties may result in some wells being unresponsive.A large noise-to-signal ratio and/or a low sampling resolution could result in low resolution measurements of the water level and any changes that occur.Monitoring well site conditions may also influence response recordings with pumping, precipitation, or large seasonal variations having an effect on water levels.Aquifers with high storage capacities and/or low transmissivity may prevent responses being recorded at the monitoring wells.This in turn may have resulted in a higher shaking threshold required for a water level change.Understanding why sites did not respond should be investigated in future work.
The characteristics of dynamic shaking are influenced by geological conditions [40].The amplification of seismic shaking is larger over sediments than bedrock [100], shown by high levels of shaking in parts of the Marlborough, Wellington, Manawatu-Wanganui, Taranaki, and Hawkes Bay regions (Figures 6 and 8).Seismic wave attenuation and velocity are also affected by the degree of saturation and the spatial distribution of fluids within the crust [101].Although we evaluated seismic and hydrogeological factors separately, we acknowledge that a nonlinear relationship exists [102] between them.

Conclusion
We quantify groundwater level changes in 433 monitoring wells across New Zealand to the 2016 M w 7.8 Kaikōura earthquake.We compare water level changes to earthquakedriven characteristics such as mean static stress changes (σ kk ), maximum peak ground velocity (PGV), and horizontal (geometric mean) peak ground acceleration (PGA).We also compare scaled water level changes (WLC/PGA) to local hydrogeological factors, depth, and site average shear wave velocity (Vs 30 ).
The Kaikōura earthquake-induced static stress changes were only significant in the near field.The amplitude and polarity of water level changes observed in monitoring wells in the near field do not generally correlate with the modelled σ kk .However, above a PGA of ~2 m/s 2 , persistent water level changes predominantly increased.This is consistent with the hypothesis of shear-induced consolidation [83].
For wells that experienced a PGA lower than ~2 m/s 2 , there was approximately an equal number of water level increases and decreases.The statistically random polarity of water level change is consistent with the hypothesis of enhanced permeability by dislodging colloids [62], as the polarity of the water level change depends on the location of the monitoring well relative to the location of the permeability change (up-or down-hydraulic head gradient) [84].
The fault's northward rupture and resulting directivity effects [5] resulted in a spatially asymmetric distribution of PGA.Water level changed persistently 60% of the time north of the epicentre, whereas south of the epicentre water level changed persistently only 48% of the time.Considering the northward directivity and that both enhanced permeability and shear-induced consolidation are primarily controlled by dynamic shaking, we find that the transition between shear-induced consolidation and enhanced permeability occurs at a PGA of ~2 m/s 2 .This threshold should be confirmed in future studies.
Hydrogeological factors depth and Vs 30 positively correlate with WLC/PGA.Regardless of the magnitude of the earthquake, monitoring wells may have a predisposition to have water level changes of certain amplitudes relative to PGA, based on the degree of confinement and the strength of the rock they are screened in.Additional work should also include changes in spring discharge and other earthquake hydrological responses to the Kaikōura earthquake, at local and catchment scales.
This immediate report examines the effect of a single earthquake on multiple hydrological sites.To provide an informative map of aquifer susceptibility to earthquakes, a multiearthquake multisite dataset, composed of individual case studies such as this one, must be developed.These datasets should collectively span significant decadal time periods and record the lack of responses as well as responses.Their analysis will distinguish the role of extrinsic (earthquakerelated) and intrinsic (local geology and hydrogeology) factors and potentially could be utilised to inform practitioners on seismic risk to aquifers and water supplies [103].

Figure 1 :Figure 2 :
Figure 1: Spatial distribution of the hydrological and seismic dataset that recorded the 2016 M w 7.8 Kaikōura earthquake.(a) Kaikōura earthquake focal mechanism [5] and surface rupture zone [13].(b) 306 seismic stations (strong motion and broadband) recorded the Kaikōura earthquake.(c) Water level data were collected from 433 monitoring wells, and rainfall and atmospheric pressure were collected from 65 climate stations, which spanned the earthquake interval.Monitoring wells in Cromwell Gorge are highlighted differently.

Figure 4 :
Figure 4: Spatial distribution of water level changes and mean static stress changes (σ kk ) induced by the Kaikōura earthquake.Red colours indicate compression; blue colours indicate tension.A transect through the surface rupture zone and the majority of monitoring wells is shown in Figure 8.

Figure 5 :Figure 6 :
Figure 5: The relationship between earthquake-induced mean static stress changes (σ kk ) and water level changes.(a) Static stress-induced compressional stress changes.High areas of compressional stress changes include (i) Wellington and Canterbury and (ii) Marlborough.(b) Static stress-induced tensional stress changes.High areas of tensional stress changes include (i) Tasman.

Figure 7 :
Figure 7: Absolute water level change amplitude as a function of (a) maximum PGV and (b) horizontal (geometric mean) PGA.Above a PGV of ~0.3 m/s and a PGA of ~2 m/s 2 , water level changes were predominantly increases.

Figure 8 :
Figure 8: A cross section of the Kaikōura earthquake-induced (a) groundwater level changes, (b) maximum PGV and horizontal (geometric mean) PGA, and (c) mean static stress changes (σ kk ).Locations of interest are (i) Cromwell Gorge monitoring wells, (ii) earthquake epicentre, (iii) Marlborough region, (iv) Cook Strait, (v) Wellington, (vi) Manawatu-Wanganui, and (vii) Hawkes Bay regions.The cross section goes through the Kaikōura surface rupture zone and the majority of monitoring wells.

Figure 9 :
Figure 9: WLC/PGA as a function of monitoring well depth.Coloured according to monitoring well type (water-table, nonflowing artesian, and flowing artesian).The deeper the monitoring well, the more sensitive the well was to water level change.All types of water level response (increase, decrease, transient, and no change) were observed over the entire depth range.Artesian wells commonly show water level increases at high WLC/PGA.

Figure 10 :
Figure10: WLC/PGA as a function of site classification[78].Site classifications are converted from Vs 30[77].Monitoring well depths shallower than 30 m have larger opaque symbols as Vs 30 is representative.Monitoring well depths deeper than 30 m have smaller transparent symbols as Vs 30 is less representative.