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 ~2 m/s2, 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 ~2 m/s2, 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.
Central New Zealand has experienced several moment magnitude (Mw) 7 or larger earthquake ruptures in the last 200 years, the largest of which were the 1848 Mw 7.4−7.7 Awatere earthquake [1, 2], the 1855 Mw 8.1+ Wairarapa earthquake , and the 1888 Mw 7−7.3 Hope earthquake . The most recent Mw 7+ event in the region was the Mw 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 , following four decades of relative quiescence [6, 7]. The Kaikōura earthquake followed the 2010-2011 Canterbury earthquake sequence, which included the devastating Mw 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 . The rupture took over two minutes , and ground acceleration exceeded ~1 g in the near-source region .
Earthquakes are known to induce groundwater level changes with varying amplitude, duration, and polarity of response (e.g., [18–20]). Earthquake-induced static and dynamic stresses (e.g., ) and local hydrogeological factors  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 earthquake-driven 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.  following the Mw 7.1 Darfield (Canterbury) earthquake, and is of comparable size to datasets compiled following the 1999 Mw 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.
2. 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 . Transgression/regression sequences created overlaps of terrestrial and marine deposits in coastal areas , resulting in large variations in vertical and horizontal permeability and storage capacity . 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 .
3. Hydrological and Seismic Data
The hydrological monitoring infrastructure in New Zealand is extensive due to the country’s heavy reliance on groundwater for agriculture . 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 , highly variable , and catchment dependent .
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 , 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–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 earthquake-induced 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 : “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), 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.50 Hz 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 . 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.
4. 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  ranging between 10–2 and 105 J/m3 (Figure 3). Although the relationship of seismic energy density, magnitude, and epicentral distance was derived from earthquakes in southern California , the term is accepted internationally and used for seismo-hydrological phenomena (e.g., Weingarten and Ge, 2013; [36, 42, 43]). Of the well water levels observed, 146 showed a persistent increase, 91 exhibited a persistent decrease, and 69 displayed a transient change. At 127 wells, there was no observed change. In the near field, approximately two-thirds of the persistent water level changes observed were increases. In the intermediate field, changes were roughly equally split between increases and decreases. The largest persistent water level increase was ~3.5 m and the largest decrease was ~3.3 m, both of which were observed in Cromwell Gorge. In simple terms, water level changes in the near field had larger amplitudes than those in the intermediate field. The time taken for water levels to reequilibrate at new postearthquake levels generally ranged from 10 minutes to two hours, with a median time of 65 minutes and the longest time of 100+ days (Cromwell Gorge).
5. The Influence of Earthquake-Driven Factors on Water Level Changes
5.1. Earthquake-Induced Static Stress Changes
Static stress changes are permanent and decay rapidly with distance () at , meaning that they are most significant in the near field (Figure 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–49]). The contractional and dilatational volumetric strain changes are inferred to increase and decrease water levels, respectively, with predictable amplitudes (e.g., ). 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 (), we used the slip distribution of Clark et al.  to calculate the internal strain field based on the formulation of Okada  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  and geodetic data. Positive indicates tensional stress changes, and negative indicates compressional stress changes.
As expected, was most significant in the near field (Figure 4). In total, 386 wells were in areas of contraction, where ranged from −0.2 to −18,800 kPa. 47 wells were in areas of dilatation, where ranged from 0.7 to 120 kPa (Figure 5). Static stress-induced compressional changes larger than 101 kPa occurred predominantly with water level increases (), and static stress-induced tensional changes between 101 and 102 kPa occurred generally with water level decreases (). Compressional changes smaller than 101 kPa occurred more frequently with decreases in water level () than increases (), although the larger amplitudes were mostly increases. At lower than ±102 kPa regardless of the sign of , some water levels did not respond. Significant (>1 m) water level changes were observed in Cromwell Gorge, where was less than −100 kPa (Figure 5).
5.2. Earthquake-Induced Dynamic Shaking
Dynamic stress changes caused by the passage of seismic waves are short-lived and decay at , meaning that they can be significant from the near to far field [44, 45]. Dynamic stress-induced deformation depends on loading rates and cycles of inertial forces  and is influenced by rupture directivity and radiation patterns .
Peak ground velocity (PGV) reflects the majority of energy in seismic ground motion  and is easily measurable. PGV can be related to the Modified Mercalli Intensity scale (e.g., ) and is a potential indicator for engineers in seismic assessment methods (e.g., ). PGV has been compared to water level changes (e.g., [35, 37]) and is related to peak dynamic stress, which has been utilised in permeability enhancement experiments (e.g., Roberts, 2005; ). Dynamic stresses are thought to control the incidence of shear-induced consolidation and dilatation [57–59], vertical (e.g., [40, 60, 61]) and horizontal (e.g., ) enhancement of permeability. Peak dynamic stress (, GPa) can be expressed  as in terms of the shear modulus ( GPa), shear wave velocity at the monitoring well (), 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., ), 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 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/s2. Above a PGV of ~0.3 m/s and a PGA of ~2 m/s2, 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/s2, 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/s2, compared to 11 persistent water level increases above ~0.3 m/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.
6. 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 .
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/s2) experienced at each monitoring well (WLC/PGA, s2). 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.
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 , as a result of the specific yield of unconfined aquifers being higher than the storativity of confined aquifers . 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 . Assuming undrained conditions , changes in confined aquifer pore pressure (, GPa), given by , are proportional to mean stress (, GPa), and the magnitude of pore pressure change is governed by Skempton’s coefficient () .
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 well was to water level change. A similar observation was recorded by Rutter et al. , 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 s2. At depths between 10 and 100 metres, WLC/PGA mostly varied from ~10–2 to 100 s2. Deeper than 100 metres, WLC/PGA largely fluctuated from ~10–1 to 101 s2. 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 100.
6.2. Shear Wave Velocity
Seismic waves attenuate differently in different rock types , and it is therefore unsurprising that hydrogeological factors partly control earthquake-induced groundwater level changes. Site effects and varying crustal structures have differing efficiencies of dynamic shaking . Saturated unconsolidated media damp high-frequency motions, while amplifying lower frequencies [25, 72, 73]. In contrast, consolidated/crystalline aquifer material amplify high-frequency shaking more than low-frequency shaking . Sensitive hydrological sites screened in well-consolidated 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 (Vs30, m/s) over depths between 0 and 30 m  at each monitoring well. Shear wave velocities range from 40 to 1040 m/s, representing very soft soil to weak rock, respectively . Vs30 is typically towards the lower end of this range, the median Vs30 being 225 m/s, representative of deep soil. At any Vs30, WLC/PGA varied over three orders of magnitude. Broadly, as Vs30 increases, WLC/PGA increases. When Vs30 ranged between 0 and 270 m/s, WLC/PGA typically ranged from ~10–3 to 100 s2. When Vs30 was close to 1000 m/s, WLC/PGA was mostly between ~10-2 and 101 s2. There is no correlation between polarity of water level change and Vs30 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 Vs30 values (Figure 10).
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 (), used in this study, can be converted into volumetric strain change (), , with Young’s modulus (, GPa) of the lithology. For Young’s modulus, an arithmetic mean of 13 GPa for sands and gravels was used . has the relationship with pore-pressure change () in a monitoring well, (e.g., ). To obtain a first-order estimate of the static strain changes and corresponding water level changes in this investigation, arithmetic means of shear modulus (, 5.6 GPa), Skempton coefficient (, 0.62), and undrained Poisson ratio (, 0.33) of sands and gravels were used , resulting in a coefficient value of ~9 GPa. Therefore, of ±100 and ±101 kPa predicts a ~10 cm and ~1 m water level change, respectively. However, most wells with of <±100 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 larger than ±100 kPa had smaller water level changes than predicted, although the polarities were generally consistent with static stress changes exceeding ±101 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–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.
7.2. Earthquake-Induced Dynamic Shaking
If dynamic stress-induced cyclic shear strains in unconsolidated deposits greatly exceed a threshold of ~10–4 , shear-induced dilatation occurs . Groundwater levels can decrease persistently as a result of the decrease in pore pressure and increase in porosity . Shear-induced dilatation may have occurred in close proximity to the Kaikōura earthquake fault rupture. However, as a result of the numerous landslides  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 . The threshold at which persistent water level increases predominantly occurred was at a horizontal (geometric mean) peak ground acceleration (PGA) of ~2 m/s2 and a maximum peak ground velocity (PGV) of 0.3 m/s. 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/s2 and had predominantly persistent water level increases.
At low thresholds of dynamic shaking, earthquake-induced flow velocities may be strong enough to dislodge colloidal particles . It has been postulated that dislodging of colloidal particles from preferential flow pathways may enhance horizontal permeability ([62, 84, 91] ). The change in water level may originate in close proximity to a local pressure source, which could be produced by liquefaction , or elevated hydraulic heads . Earthquake-induced flow velocities may be strong enough to dislodge colloidal particles . Laboratory experiments on pore unclogging [56, 94] and groundwater colour changes  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 . In this study, there was roughly an equal number of persistent water level increases and decreases at accelerations lower than PGA of ~2 m/s2, 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 with the model of enhanced permeability. Persistent water level changes may have occurred above a PGA of ~2 m/s2 as a result of enhanced permeability but may be concealed by larger amplitude changes caused by shear-induced consolidation .
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 . In this study, the median persistent response reequilibration time for monitoring wells that experienced a PGA above ~2 m/s2 () was 585 minutes, while below a PGA of ~2 m/s2 () was 75 minutes. The two populations of response times were statistically different. Shear-induced consolidation may occur on a large scale (km) , 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)  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 . In the field, seismic energy densities larger than ~101 J/m3 are above the threshold for undrained consolidation, with sensitive sites experiencing consolidation above ~10–1 J/m3 [38, 96]. For the Chi-Chi and Hengchun earthquakes, the transition from consolidation to enhanced permeability was inferred at ~101 J/m3 . 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 , near and intermediate fields , and one fault rupture length  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/s2 (Figure 7), where liquefaction occurred, and the random polarity of responses below ~2 m/s2, we find that the transition between shear-induced consolidation and enhanced permeability may occur at a PGA of ~2 m/s2. 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  and poroelastic deformation (e.g., ). In this case study, transient water level changes occurred mainly below a PGA of ~2 m/s2 with small amplitudes of water level changes, (Figure 7), and are likely to be a result of poroelastic deformation (e.g., [64, 92]).
7.3. 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 . 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/m3 . Since there is a disparity between the global dataset  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., ).
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 . Also, with depth, the transmissivity generally becomes lower and the strain sensitivity of the water head becomes larger .
WLC/PGA is positively correlated with Vs30 (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.  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 Vs30. 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 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 . The amplification of seismic shaking is larger over sediments than bedrock , 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 . Although we evaluated seismic and hydrogeological factors separately, we acknowledge that a nonlinear relationship exists  between them.
We quantify groundwater level changes in 433 monitoring wells across New Zealand to the 2016 Mw 7.8 Kaikōura earthquake. We compare water level changes to earthquake-driven characteristics such as mean static stress changes (), 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 (Vs30).
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 . However, above a PGA of ~2 m/s2, persistent water level changes predominantly increased. This is consistent with the hypothesis of shear-induced consolidation .
For wells that experienced a PGA lower than ~2 m/s2, 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 , 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) .
The fault’s northward rupture and resulting directivity effects  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/s2. This threshold should be confirmed in future studies.
Hydrogeological factors depth and Vs30 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 (earthquake-related) and intrinsic (local geology and hydrogeology) factors and potentially could be utilised to inform practitioners on seismic risk to aquifers and water supplies .
The data supporting the results reported can be found in Supplementary Materials.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
We would like to thank many organisations and people for providing data for this study: Northland Regional Council (Sandrine Le Gars, Alan Bee, and Susie Osbaliston); Auckland Council (Nicholas Holwerda); Bay of Plenty Regional Council (Diane Harvey, Brent Hutchby); Gisborne District Council (Matthew McGill-Brown, Peter Hancock, and Murry Cave); Hawkes Bay Regional Council (Simon Harper); Waikato Regional Council (John Hughey); Taranaki Regional Council (Jane Harvey, Fiona Jansma, and Regan Phipps); Horizons Regional Council (Stephen Collins, Brent Watson); Greater Wellington Regional Council (Sheree Tidswell, Doug Mzila, and Mike Thompson); Marlborough District Council (Peter Davidson); Tasman District Council (Joseph Thomas, Monique Harvey); Environment Canterbury (Shaun Thomsen); Otago Regional Council (Andrew Egan, Nineva Vaitupu); Environment Southland (Michael Killick); Watercare (Andrew Lester); GNS Science (Abigail Lovett, Grant O’Brien); Contact Energy (Neil Whitford); Golder Associates (Eric van Nieuwkerk); Aqualinc Research Ltd. (Ross Hector); Wairakei Estate (Nic Conland); Manfeild (Steve Easthope); Victoria University of Wellington (Rupert Sutherland); and the National Institute of Water and Atmospheric Research. We also wish to thank colleagues Mai-Linh Doan and Guðjón Eggertsson for their discussion and helpful comments. We acknowledge the New Zealand GeoNet project and its sponsors EQC, GNS Science, and LINZ for providing data used in this study. This study was funded under the Royal Society of New Zealand Marsden Fund (2012-GNS-003).
Supplementary data for each monitoring well include the earthquake-induced water level change measurements and the corresponding seismic and hydrogeological variable values. The supplementary figure shows groundwater level relative to ground against depth of monitoring well with symbols representing response type (increase, decrease). (Supplementary Materials)
D. W. Rodgers and T. A. Little, “World’s largest coseismic strike-slip offset: the 1855 rupture of the Wairarapa Fault, New Zealand, and implications for displacement/ length scaling of continental earthquakes,” Journal of Geophysical Research, vol. 111, no. B12, 2006.View at: Publisher Site | Google Scholar
G. Downes and D. Dowrick, Atlas of Isoseismal Maps of New Zealand Earthquakes, 1843-2003, Second Ed. Technical Report, Institute of Geological & Nuclear Sciences, Lower Hutt, New Zealand, 2015.
I. J. Hamling, E. D'Anastasio, L. M. Wallace et al., “Crustal deformation and stress transfer during a propagating earthquake sequence: the 2013 Cook Strait sequence, central New Zealand,” Journal of Geophysical Research: Solid Earth, vol. 119, no. 7, pp. 6080–6092, 2014.View at: Publisher Site | Google Scholar
C. Holden, A. Kaiser, R. Van Dissen, and R. Jury, “Sources, ground motion and structural response characteristics in wellington of the 2013 Cook Strait earthquakes,” Bulletin of the New Zealand Society for Earthquake Engineering, vol. 44, no. 4, pp. 188–195, 2013.View at: Google Scholar
C. Holden, Y. Kaneko, E. D'Anastasio, R. Benites, B. Fry, and I. J. Hamling, “The 2016 Kaikōura earthquake revealed by kinematic source inversion and seismic wavefield simulations: slow rupture propagation on a geometrically complex crustal fault network,” Geophysical Research Letters, vol. 44, no. 22, pp. 11,320–11,328, 2017.View at: Publisher Site | Google Scholar
M. W. Stirling, N. J. Litchfield, P. Villamor et al., “The Mw7.8 2016 Kaikoura earthquake: surface fault rupture and seismic hazard context,” Bulletin of the New Zealand Society for Earthquake Engineering, vol. 50, no. 2, pp. 73–84, 2017.View at: Google Scholar
B. A. Bradley, H. N. T. Razafindrakoto, and M. A. Nazer, “Strong ground motion observations of engineering interest from the 14 November 2016 Mw 7.8 Kaikoura, New Zealand earthquake,” Bulletin of the New Zealand Society for Earthquake Engineering, vol. 50, no. 2, pp. 85–93, 2017.View at: Google Scholar
M. Petitta, L. Mastrorillo, E. Preziosi et al., “Water-table and discharge changes associated with the 2016-2017 seismic sequence in central Italy: hydrogeological data and a conceptual model for fractured carbonate aquifers,” Hydrogeology Journal, vol. 26, no. 4, pp. 1009–1026, 2018.View at: Publisher Site | Google Scholar
G. Tranfaglia, E. Esposito, S. Porfido, and R. Pece, “The 23 July 1930 earthquake (Ms= 6.7) in the Southern Apennines (Italy): geological and hydrological effects,” Bollettino Geofisico, vol. 1, no. 1, pp. 63–86, 2011.View at: Google Scholar
A. Besedina, E. Vinogradov, E. Gorbunova, and I. Svintsov, Chilean Earthquakes: Aquifer Responses at the Russian Platform, In the Chile-2015 (Illapel) Earthquake and Tsunami (Pp. 133-144), Birkhäuser, Cham, 2017.
P. A. White, Groundwater Resources in New Zealand, Groundwaters of New Zealand, 2001.
A. Bal, “Valley fills and coastal cliffs buried beneath an alluvial plain: evidence from variation of permeabilities in gravel aquifers, Canterbury Plains, New Zealand,” Journal of Hydrology, New Zealand, vol. 35, no. 1, pp. 1–27, 1996.View at: Google Scholar
G. A. O'Brien, S. C. Cox, and J. Townend, “Spatially and temporally systematic hydrologic changes within large geoengineered landslides, Cromwell Gorge, New Zealand, induced by multiple regional earthquakes,” Journal of Geophysical Research: Solid Earth, vol. 121, no. 12, pp. 8750–8773, 2016.View at: Publisher Site | Google Scholar
C. Y. Wang and M. Manga, Earthquakes and Water, Springer, 2010.
D. Ebdon, Statistics in Geography, Blackwell, 1985.
G. B. Cua, Creating the Virtual Seismologist: Developments in Ground Motion Characterization and Seismic Early Warning, [Ph.D Dissertation], California Institute of Technology, 2004.
S. C. Cox, C. D. Menzies, R. Sutherland, P. H. Denys, C. Chamberlain, and D. A. H. Teagle, “Changes in hot spring temperature and hydrogeology of the Alpine Fault hanging wall, New Zealand, induced by distal South Island earthquakes,” Geofluids, vol. 15, no. 1-2, 239 pages, 2015.View at: Publisher Site | Google Scholar
T. Lay and T. C. Wallace, Modern Global Seismology, Academic Press, San Diego, 1995.
E. Quilty and E. A. Roeloffs, “Water-level changes in response to the December 20, 1994, M4.7 earthquake near Parkfield, California,” Bulletin of the Seismological Society of America, vol. 87, no. 2, pp. 310–317, 1997.View at: Google Scholar
Y. Okada, “Internal deformation due to shear and tensile faults in a half-space,” Bulletin of the Seismological Society of America, vol. 82, no. 2, pp. 1018–1040, 1992.View at: Google Scholar
M. C. Gerstenberger, C. B. Worden, and D. J. Wald, “A probabilistic relationship between ground shaking parameters and MMI based on felt report data,” in Conference Proceedings, New Zealand Society for Earthquake Engineering, Palmerston North: Performance by design, can we predict it? Palmerston North, New Zealand, 2007.View at: Google Scholar
R. Dobry, R. S. Ladd, F. Y. Yokel, R. M. Chung, and D. Powell, “Prediction of pore water pressure buildup and liquefaction of sands during earthquakes by the cyclic strain method,” in National Bureau of Standards Building Science Series, 138, p. 150, National Bureau of Standards and Technology, Gaithersburg, MD, USA, 1982.View at: Google Scholar
M. P. Luong, “Stress-strain aspects of cohesionless soils under cyclic and transient loading,” in Proc. Intern. Symp. Soils under Cyclic and Transient Loading, vol. 1, pp. 315–324, A. Balkema, Rotterdam, Netherlands, 1980.View at: Google Scholar
J. Jaeger and N. Cook, Fundamentals of Rock Mechanics, Chapman and Hall, London, 3rd edition, 1979.
A. E. H. Love, Elasticity, Dover, New York, 1927.
B. A. Bradley and M. Hughes, Conditional Peak Ground Accelerations in the Canterbury Earthquakes for Conventional Liquefaction Assessment, University of Canterbury, Christchurch, New Zealand, 2012.
New Zealand Geotechnical Society, Geotechnical Earthquake Engineering Practice: Module 1 – Guideline for the Identification, Assessment and Mitigation of Liquefaction Hazards, p. 34, 2010.
R. A. Freeze and J. A. Cherry, Groundwater, Prentice-Hall, Inc, Englewood Cliffs, NJ, USA, 1979.
N. A. Horspool, J. Chadwick, J. Ristau, J. Salichon, and M. C. Gerstenberger, “ShakeMapNZ: informing post-event decision making,” in 2015 NZSEE Conference, pp. 370–376, Rotorua, New Zealand, 2015.View at: Google Scholar
U. Destegul, G. Dellow, and D. W. Heron, “A ground shaking amplification map for New Zealand,” Bulletin of the New Zealand Society for Earthquake Engineering, vol. 42, no. 2, p. 122, 2008.View at: Google Scholar
S. Dellow, C. Massey, S. Cox et al., “Landslides caused by the Mw7.8 Kaikoura earthquake and the immediate response,” Bulletin of the New Zealand Society for Earthquake Engineering, vol. 50, no. 2, pp. 106–116, 2017.View at: Google Scholar
A. J. Davies, V. Sadashiva, M. Aghababaei et al., “Transport infrastructure performance and management in the South Island of New Zealand, during the first 100 days following the 2016 Mw 7.8 “Kaikoura” earthquake,” Bulletin of the New Zealand Society for Earthquake Engineering, vol. 50, no. 2, pp. 271–299, 2017.View at: Google Scholar
M. E. Stringer, S. Bastin, C. R. McGann et al., “Geotechnical aspects of the 2016 Kaikoura earthquake on the South Island of New Zealand,” Bulletin of the New Zealand Society for Earthquake Engineering, vol. 50, no. 2, pp. 117–141, 2017.View at: Google Scholar
M. Cubrinovski, J. D. Bray, C. De La Torre et al., “Liquefaction effects and associated damages observed at the Wellington centreport from the 2016 Kaikoura earthquake,” Bulletin of the New Zealand Society for Earthquake Engineering, vol. 50, no. 2, pp. 152–173, 2017.View at: Google Scholar
R. P. Orense, Y. Mirjafari, S. Asadi et al., “Ground performance in Wellington waterfront area following the 2016 Kaikoura earthquake,” Bulletin of the New Zealand Society for Earthquake Engineering, vol. 50, no. 2, pp. 142–151, 2017.View at: Google Scholar
N. Matsumoto, G. Kitagawa, and E. A. Roeloffs, “Hydrological response to earthquakes in the Haibara well, central Japan - I. Groundwater level changes revealed using state space decomposition of atmospheric pressure, rainfall and tidal responses,” Geophysical Journal International, vol. 155, no. 3, pp. 885–898, 2003.View at: Publisher Site | Google Scholar
J. C. Prior and P. J. Lohmann, Iowa’s Groundwater Basics, Iowa Department of Natural Resources, 2003.
K. Hazirbaba and E. M. Rathje, “A comparison between in situ and laboratory measurements of pore water pressure generation,” in 13th World Conference on Earthquake Engineering, Vancouver, 2004.View at: Google Scholar
D. L. Wells and K. J. Coppersmith, “New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement,” Bulletin of the Seismological Society of America, vol. 84, no. 4, pp. 974–1002, 1994.View at: Google Scholar
J. P. Eaton and K. J. Takasaki, “Seismological interpretation of earthquake induced water-level fluctuations in wells,” Bulletin of the Seismological Society of America, vol. 49, no. 3, pp. 227–245, 1959.View at: Google Scholar
H. T. Stearns, “Record of earthquake made by automatic recorders on wells in California,” Bulletin of the Seismological Society of America, vol. 18, no. 1, pp. 9–15, 1928.View at: Google Scholar
P. Y. Bard and J. Riepl-Thomas, “Wave propagation in complex geological structures and their effects on strong ground motion,” in Wave Motion in Earthquake Engineering, E. Kausel and G. D. Manolis, Eds., pp. 37–95, International Series Advances in Earthquake Engineering, WIT, Boston, MA, USA, 1999.View at: Google Scholar
I. A. Beresnev and K. L. Wen, “Nonlinear soil response—a reality?” Bulletin of the Seismological Society of America, vol. 86, no. 6, pp. 1964–1978, 1996.View at: Google Scholar
K. C. Weaver, Earthquake-Induced Hydrogeological Changes in New Zealand, [Ph.D Thesis], Victoria University of Wellington, 2018.