Mapping the pre and post seismic crustal stress heterogeneity analogous to 2016, Mw 7.8, Kaikoura quake, New Zealand

Heterogeneity of pre and post seismic stress states associated to any earthquake play a primary role in understanding the earthquake mechanism and hazard assessment of a seismically dynamic region. The M w 7.8, November 14, 2016 Kaikoura, New Zealand earthquake offer an unprecedented possibility to observe the heterogeneity in stress eld over a very complex fault system wherein subduction zone converges with the strike slip faults system. Here we report the pre and post seismic stress eld asperity rst time in terms of spatial and temporal variations of b-values associated to the Kaikoura main-shock. Pre seismic spatial disparity of b-value indicates the existence of two prominent low b-value clusters, one towards southwest closer to the epicenter and other to the north of the rupture zone. During co seismic period, owing to the stress release near the epicentral area, the pattern of prominent low b-value pattern has become negligible in the post seismic period. However, the pattern of low b-value in the north of the rupture zone remains similar in the post seismic period indicates the unreleased strain energy in the province. The temporal evaluation of the earthquakes frequency magnitude distributions over a period of two decades also showed an analogous pattern that the b-values were decreased considerably before the large earthquakes in the expanse, which could spawn a larger future earthquakes in the vicinity. Continuous black lines indicate the active fault systems and subduction zones of the region. The white arrow represents the co-seismic rupture propagation direction towards northeast from the epicenter.


Introduction
On November 14, 2016 at 12:03 a.m. local time, an effective earthquake of magnitude M w 7.8 rattled the Kaikoura vicinity within the northeastern coastal area of south Island, New Zealand. The earthquake struck in the tectonically more complex region, in which a transition of Hikurangi subduction zone converge with a strike-slip dominated Alpine Fault via Marlborough fault systems within the continental area of northeastern province of the South Island ( Fig. 1) (Kaiser et al., 2017). The event triggered at a shallow depth of 15 km, propagated the rupture northeastward about 170 km across the New Zealand plate boundary zone from south of eastern Hope Fault to Cape Campbell and Cook Strait (Duputel & Rivera, 2017;Hamling et al., 2017;Nicol et al., 2018). The focal mechanism indicated a non-double couple source, with signi cant right-lateral strike slip and thrust motion. Within three days, more than 2000 aftershocks took place associated to the main-shock in which four where M w > 6. Consequently, the event gained the eye of geoscienti c netwrok for its complex surface rupture and mechanism across ~ 12 distinct active fault segments in two speci c tectonic domains, mainly the strike-slip Marlborough Fault System (MFS) and the transpressional North Canterbury Domain (NCD) (Hamling et al., 2017;Kaiser et al., 2017).
The number of earthquakes as a function of their sizes had been considered to follow a self-comparable power-law distribution. One such power-law relation is the frequency-magnitude relation described by using the b-value proposed by Gutenberg and Richter in California (1944) and by Ishimoto and Iida (1939) in Japan. The b-value in the Gutenberg-Richter relation, as it is commonly called, is usually close to 1 over long time periods and large expanses (Lay and Wallace, 1995). However, statistically signi cant variations have been observed for shorter time windows and limited geographical regions. In addition to the stress state of the different parts of medium, the b-value re ects the connection among large and small earthquakes. According to Urbancic (1992), the size distribution of earthquakes quanti ed by bvalue is possibly a proxy to differentiate stress conditions and could consequently act as a crude stressmeter wherever seismicity is observed. Thus, the b-value is often interpreted as an indicator of applied shear stress and material heterogeneity. Several studies (Wyss et al., 2000;Schorlemmer et al., 2004;Nuannin et al., 2005;Wyss and Stefansson, 2006;Ghosh et al., 2008;Nuannin et al.,2012) have proven that spatial and temporal b-value anomalies precede the prevalence of large earthquakes. For the reason that Kaikoura earthquake of M w 7.8 was a major event in the region with a complex tectonic mechanism, the b-value anomalies concomitant with this event is of great interest in terms of earthquake investigation. These anomalies, if any present, will shed light on the stress distribution of the region.
Numerous studies (Scholtz, 1968;Wyss, 1973;Wu et al,. 2008b;Nanjo et al., 2012) suggest a decrease in b-value with durations of months and years interpreted as boom in stress prior to the occurrence of a main event. This could mean that low b-values can be expected near the hypocenters of future earthquakes. In this study, we attempt to probe the variation of Gutenberg-Richter b-value, both spatial and temporal, and explore the disparity of the pre and post seismic crustal stress conditions associated to 2016 Mw 7.8 Kaikoura earthquake, northeastern south Island of New Zealand (Fig. 1).

Tectonic Settings Of Central New Zealand
The country of New Zealand straddles the boundary between the active obliquely convergent Paci c-Australian plates which converge at rates of 39-48 mm/yr (Beavan et al., 2002). The opposite-dipping Hikurangi subduction systems linked by the Alpine Fault transform through MFS provides a su ciently large seismicity across New Zealand. In a comparatively short, documented period, the numerous active faults associated with these plate boundaries have produced frequent large earthquakes in New Zealand.
The relative plate motion changes from slightly oblique subduction of the Paci c Plate beneath the Australian Plate at a rate close to 50 mm/yr in the northeast of north Island to transpression and very oblique continental collision at a rate < 40 mm/yr in the south Island (DeMets et al., 2010).
The M w 7.8, 2016 Kaikoura earthquake occurred within the area of transition between the southern Hikurangi Subduction Zone and the oblique continental collision along the Alpine Fault (Fig. 1). The MFS and NCD collectively mark this transition. The MFS includes four major right-lateral strike-slip faults from north to south, the Wairau, Awatere, Clarence and Hope faults. The faults within this system head southwards and converge to form the Alpine fault. The Alpine fault, a 650-km-long dextral strike-slip fault is a main tectonic feature of the south Island. Several dextral faults cut through the upper South Island and much of the north Island. The focal mechanism solution (U.S. Geological Survey, 2016) indicates that the preferred fault plane of the strikes northeast roughly parallel to the coastline with a down-dip toward northwest direction. The event occurred along the multi-segmented faults of around 12 faults in a combination of different orientations and slip mechanisms. The rupture of the event has propagated from southwest to northeast along with strike-slip extending over the rupture length of ~170 km from the epicenter (Duputel & Rivera, 2017;Hamling et al., 2017;Nicol et al., 2018;Xu et al., 2018). The aftershock zone comprises of the epicenter near the Humps Fault zone and progressing northeastward through the Hundalee Fault, the Hope fault, the Jordan thrust, the Papatea Fault and the Kekerengu fault before stopping at the Cook Strait. The region has seen experienced several large M > 7 earthquakes in the past like the 1848 Marlborough earthquake, 1855 Wairarapa earthquake, 1888 Amuri earthquake 1929 Arthurs pass earthquake etc. Several earthquakes of magnitudes between 6 and 7 has been recorded in the region too.

Data And B-value Mapping
The earthquake size distribution in the Earth's crust commonly follows the Gutenberg-Richter power law (Gutenberg & Richter, 1944). The Gutenberg-Richter law gives a fundamental statistical description of seismicity as: log 10 N = a -bM, where a and b are constants, M is the magnitude, and N is the number of earthquakes that occur in a speci c time window with magnitude ≥ M. According to this law, a lower bvalue means that larger earthquakes will dominate over smaller earthquakes and a higher b-value means that there is dominance of small earthquakes over large earthquakes.
The completeness of the earthquake catalog used is important for accurate b-value estimations (Nuannin, et al., 2012). This is achieved by the estimation of magnitude of completeness (M c ) which is the smallest magnitude at which the catalog is thought to have included all seismic events. An underestimate of M c can result in too low b-value while a larger M c value increases the uncertainty in bvalue estimation. Among different techniques for M c calculation, M c for this catalog is calculated using maximum curvature method (Wiemer&Wyss,2000). To remove the dependent events like foreshocks and aftershocks from the catalog, declustering is done. The Reasenberg(1985) algorithm is used for this purpose.
The b-value is determined by maximum-likelihood method (Aki,1965). To map the spatial variation of b, earthquake epicenters are projected onto a plane. The b-value is estimated at every node of a plane grid, using the N (N = constant) nearest earthquakes located within a distance, R, from the node. In this method, we can ensure a constant sample size for each calculation. The variation of b-value is identi ed by creating dense spatial grids and sampling overlapping circular volumes in space and plotting them on a map. The grid spacing is chosen accordingly so that there is signi cant overlapping of volumes and provide a natural smoothening of results.
Temporal b-value variations are calculated by maximum curvature method using sliding time windows containing a constant number of events. Windows with constant number of events is better than that of constant time width so that a change in sample size does not affect the calculations.
New Zealand's GNS Science has been monitoring seismic networks all over the region and it provides earthquake data. The earthquake catalog for this study was taken from GeoNet service of GNS Science for the time-period 1st January 2000 to 1st July 2020. It consisted of 17231 events with magnitudes between 3 and 7.8 and depth from 0 to 350 km. Catalog consisted of magnitude types M w , m b and M L . To make the catalog homogeneous, m b and M L were converted to M w .
The catalog has been divided into two sub-catalogs: rst one with events before the 2016 event beginning from January 1, 2000 (Fig. 2a) and second one with earthquakes after the event until July 1, 2020 (Fig. 2b). The number of events used for the pre seismic stress eld analysis is much higher than the data used for the post seismic analysis. However the number of events we have considered for both the cases are su cient for the estimation of b-value calculations (Zhao and Wu, 2008;Meng et al., 2018).
The calculations are carried out using ZMAP software package (Wiemer, 2001).

Results And Discussion
It is globally established that b-value from the Gutenberg-Richter relationship decreases with increasing differential stress level (Scholz, 1968;Scholz, 2015). To understand the stress con guration of the study area, we estimate the b-value variations before and after the 2016 Kaikoura earthquake. We apply the maximum-likelihood method, which has been proved to be a robust and unbiased method in most cases (Wiemer 2001), to compute the b-value to events in the study area. The selected earthquake catalog is closely examined to analyze the seismicity of the study area. We used magnitude values from 3 to 8 to obtain the seismicity map of the entire data catalogue for the study area. Histograms of depth and magnitude are plotted for the entire catalog ( Fig. 3a & 3b). The depth of the events considered vary from 0 to 350 km. It may be noted that maximum number earthquakes are located below 20 km depth. However, the number of earthquake events maximum up to a depth level of ~ 250 km indicates the deeper level of subduction zone seismic activities in the plate boundary of the study region. The histogram of magnitude distribution earthquake as a function of time shows that majority of the earthquake magnitudes are ranging between M w 3-4 with large magnitudes are fewer in number.
The cumulative rate of events plotted with time in Fig. 4. The cumulative number shows earthquake events as function of time. It is noticed that the rate of earthquakes was not constant through time. The earthquakes with highest magnitudes are mostly occurred in the year 2010 (M w 7.2) and 2016 (M w 7.8) where sharp increases in the pattern of cumulative rates are present. The cumulative frequency of earthquakes as a function of magnitude is shown in Fig. 5. The straight line represents the maximum likelihood estimate of a and b values based on the events above magnitude of completeness M c is determined. The frequency-magnitude distribution (FMD) for the complete catalog showed magnitude of completeness M c as 3.1, b-value as 0.87 and a-value as 6.837. All the events with magnitude less than M c were removed from the catalog. After removing events below M c and declustering, the resulting catalog consisted of 8239 events. This nal catalog is used throughout the seismicity analysis in this study.

Spatial variation of b-value
To examine the spatial distribution of the pre and post seismic b-value, the study area is divided into a dense spatial grid of 0.05 0 x 0.05 0 . b-value is calculated for circular epicentral areas centered at grid nodes where the radii of circular regions vary along the grid to cover a xed number of events. In this case, a circular region containing a minimum of 150 nearby events are required to calculate b-value in each node. The b-values are calculated by using the Maximum Likelihood method (Aki, 1965). The calculated pre and post seismic b-values are color coded and plotted as maps in Figs. 6a and 7a. The b-value map in Fig. 6a shows the pre seismic variation of b-values for a period from the year 2000 to the day before Mw 7.8 Kaikoura earthquake main-shock. The sub-catalog used for this b-value estimation consisted of 6809 events (Fig. 2a). The event epicenter and rupture zone associated to the Kaikoura event are shown in the map. It is evident that two b-value zones are clustered southwest (Zone-I) and northern (Zone-II) portions of the Kaikoura earthquake rupture zone, suggesting the high stress buildup prior to the main-shock closer to epicenter and northern end of the rupture respectively. The b-value estimated using FMD on the nodes placed over this southwest and northern stress elds are 0.78 ± 0.07 and 0.73 ± 0.07 respectively (Figs. 6b and 6c). Figure 7a shows the estimated b-value map for the post seismic period of the Kaikoura event (i.e from 15th of November 2016 to July 2020). Although the number of events after the 2016 event is comparatively lower (1427 events) than that of the pre seismic data, they are su cient to compute the b-values (Zhao and Wu, 2008;Meng et al., 2018). During the post seismic period, the bvalue of the epicentral node is found to have increased to 1.01 ± 0.08 (Fig. 7b). However, in the northern part of the rupture zone where the second patch of high stress eld demarcated during the pre-seismic period, the b-value remain same as 0.73 ± 0.05 (Fig. 7c) during the post seismic period as well.
It is established that, before an earthquake low b-value correlates with high stress in the unruptured area. Vice versa, after an earthquake b-value increases owing to stress releasing in the rupture zone . Our epicentral mapping of relatively low b-value prior and increase of b-value later to the Kaikoura earthquake is in good agreement with the studies reported by Wyss et al.  (Wiemer and Wyss, 1997;Sreejith et al., 2018). Figure 6a shows the southwest cluster of high stress eld (Zone-I) near the epicentral area located over the Marlborough fault systems in the northeastern part of the south Island (Fig. 1). It may be noted that, the triggering of Kaikoura earthquake was very closely associated to one of the extremely complex mechanism ever recorded, which involved ~ 12 separate fault segments extending over ~ 170 km from the epicenter in north Canterbury domain to Cook Strait (Kaiser et al., 2017). Investigation using the eld observations in conjunction with interferometric synthetic aperture radar, Global Positioning System, and seismology data reveal the propagation of the rupture from south Hope fault to north Cook Strait region (Hamiling et al., 2017). It is also reported that the northward propagation of the rupture along these numerous faults was the result of static stress stored over a period of time during the interseismic period due to the right-lateral oblique-slip movements in north Canterbury and right-lateral strike-slip motions in the Marlborough fault system (Jiang et al., 2018). Thus our pre seismic analysis result supports the ndings, which state that large areas of high stress accumulation were formed when all of these oblique and right lateral strike slip motions on crustal faults were employed to stress the shallow fault interface source of north Canterbury domain and Marlborough fault systems.
However after the earthquake the high pre seismic stress eld demarcated is has changed to comparatively a stress free zone due to the stress release during the co seismic rupturing (Fig. 7a). The co seismic deformation taken place not only along the shallow crustal faults of the Marlborough fault systems but also to some extent at the Hikurangi subduction interface, which slowly released the energy (Jiang et al., 2018). However in both cases, there is a clear existence of high stress eld (low b-value) in the Zone-II located over the north end of the rupture zone at Cook Strait closer to the Wellington area in southwestern part of the north Island. Jiang et al. (2018) reported from GPS observations that no effective co seismic and post seismic deformations associated with the Kaikoura event encouraged to release the accumulated stress in Cook Strait and neighboring Wellington area.
The Cook Strait has been considered a potential region for frequent earthquakes by many researchers (Wallace et al., 2009;Pondard and Barnes 2010). The tectonics of Cook Strait is mainly controlled by the transition of subduction in the north island to transpression in the northern South Island through numerous fault lines. Using GPS, Seismological and Geological data Wallace et al. (2009) reported that, Cook Strait province is having a potential to trigger an earthquake of ~ M w 8.5. However, during 2013 M w 6.6 Cook Strait earthquake and allied events sequence, released on only a part of the reported energy (Pondard and Barnes 2010). Thus a good correlation of strain rate obtained from the GPS observations (Jiang et al., 20018) and b-value estimated in the present study clearly indicate an ongoing strain buildup in the Cook Strait and Wellington regions which may trigger a larger earthquake in future.

Temporal variation of b-value
Several laboratory rock-fracture experiments (Mogi, 1963;Scholz, 1968;Lei, 2003) indicate a systematic decrease in b-value approaching the time of entire fracture. Investigation by Henderson et al., (1994) in the southern Californian region also supports this result, that seismicity preceding large main-shocks can be associated with low b-values. In the light of these results, we also evaluated the temporal variation of b-values (Fig. 8) in terms of earthquakes FMD in the study region for the past two decades. The b-values are calculated in successive time windows containing a constant number of events N. The window is moved in time by 10% increments of time. A moving window of 250 events with a minimum of 50 events and an overlap of 10% increments of events is used.
The events considered span for the period January 1, 2000 to July 1, 2020 with M > 3. In this period, two large earthquakes of magnitude M > 7 has occurred in the region. The b-value in the region varies between 0.7 and 1.2 and this range of values are consistent with the b-values obtained in the spatial analysis. We can see that before each M > 7 event that occurred in the region, the trend of the b-value decreases. Also the b-value of the region has remained below 1 for most of the major events in the past two decade suggesting that the region was under considerable stress building during this period. The lowest b-value (~ 0.71) is observed just before the 2016 Kaikoura earthquake for the entire region. It is known that before any major earthquakes, the b-value is considerably below 1. In a usual case, it is observed that the b-value decreases before a main event and then slightly increases after the event . It may be noted that, before the M w 7.8 Kaikoura event, the largest earthquake was M w 7.2 Canterbury earthquake on 3rd September 2010 (Atzori et al., 2012). Figure 8 portrays that, there was a slight decrease in b-value below one prior to Canterbury main event (Shcherbakov et al., 2012). Conversely there was a slight increase in b-value after the 2010 event. However, the average value remains well below 1 and keeps decreasing to mark the lowest b-value (~ 0.71) just before the major event of M w 7.8 in 2016. In the same mode, there is also a slight increase in the b-value after the Kaikoura earthquake. However, it is similarly noticeable that the b-value in the region has not risen above 1 after the Kaikoura event. Thus, this could be a sign of an impending earthquake towards Cook Strait province. This is consistent with the low b-value determined in the spatial variation map of the region after the Kaikoura main-shock (Fig. 7). Thus the change in b-value can be used as a stress-meter for an impending earthquake in future (Zang et al., 1998;Lei, 2003;.

Conclusions
The slope of the relation among frequency-magnitude distributions (FMD) is theoretically 1, however it often indicates deviations from 1 relies upon the heterogeneity of the area. In order to compare the pre and post seismic stress variations of in the context of the complex rupturing and mechanism of 2016 M w 7.8 Kaikoura earthquake, spatial and temporal variations of the earthquake b-values earlier and later to the main-shock are analyzed. The FMD for the entire catalog con rmed magnitude of completeness M c as 3.1, b-value as 0.87 and a-value as 6.837. The analysis of the pre seismic data envisioned a low bvalue with the vicinity of the epicenter at the outer edge of the cluster advocate a high stress accumulation prior to the main event. Conversely, the presence of a subsequent low b-value cluster identi ed to the north of Kaikoura rupture zone during pre and post seismic periods portrays an unreleased strain energy after the Kaikoura quake rupturing. The assessment of the temporal variations of b-values for almost two decades strengthened the consistency of the spatial b-value patterns before and after the major earthquake events in the study region propose that the change in b-value can be used as a stress-meter for an impending earthquake in destiny.

Declaration of competing interest
The authors declare that they have no known competing nancial interests or personal relationships that could have appeared to in uence the work reported in this paper.