Piezomagnetic Anomalies Associated with the 2021 M W 7.3 Maduo (China) Earthquake

: Stress changes due to earthquake rupture can disturb geomagnetic ﬁelds signiﬁcantly. In order to investigate the impact of the 2021 M W 7.3 Maduo earthquake on geomagnetic ﬁelds, a piezomagnetic model is constructed based on the coseismic slip to calculate the static coseismic piezomagnetic anomalies (PMs). The PMs are considerable in near-ﬁeld. However, the PMs are negligible in regions tens of kilometers from the fault rupture. The PMs of our model are consistent with those of other strike-slip earthquakes, indicating that our piezomagnetic model is reasonable. The east component of observed coseismic geomagnetic changes and calculated PMs on a geomagnetic repeat station located about 6 km from fault trace are +4.8 ± 2.2 nanotesla and +4.3 nanotesla, respectively. It seems that the piezomagnetic model can explain the observed data. The PMs are up to 10 nanotesla in the near-ﬁeld with the initial magnetization of 3 A/m and stress sensitivity of 2 × 10 − 3 MPa − 1 . Consequently, considerable coseismic geomagnetic changes that are above error could be observed along the fault, especially at locations with geometrical complexities.


Introduction
Seismomagnetic anomalies refer to magnetic disturbances associated with seismic activities [1,2]. Seismomagnetic anomalies have been observed in recent decades by highprecision magnetometers in America [3][4][5][6], Japan [7][8][9], China [10,11], etc. Seismomagnetic anomalies may provide valuable insights into stress accumulation in the crust [12] and the mechanisms of earthquakes [10]. Many seismomagnetic anomalies are coseismic changes in the geomagnetic field [7,13,14]. Others are preseismic anomalies such as anomalies observed almost simultaneously with the minima of the variability of the order parameter of seismicity [15] two to three months before major earthquakes [16,17]. One mechanism of coseismic geomagnetic changes is the piezomagnetic effect caused by the static stress changes due to earthquake rupture [18]. Many observed seismomagnetic anomalies have been explained by the piezomagnetic effect [19,20]. The piezomagnetic effect can also be used to explain the secular variation of geomagnetic fields in tectonic active zones, such as subduction zones [21]. The piezomagnetic anomalies (PMs) associated with strong earthquakes in low strain rate regions are more important because coseismic PMs are much larger than those resulting from aseismic tectonic loading [22][23][24]. Thus, coseismic PMs must be investigated and excluded before examining geomagnetic field variations due to aseismic tectonic loadings. The piezomagnetic effect can also contribute to understanding the mechanisms of earthquakes as well as their rupture processes, since PMs are significantly dependent on fault geometry, slip, etc. [7,25]. New research into PMs has also been involved in the early warning of earthquakes, because PMs may be observed earlier than seismic waves [7,26].
The Bayan Har block, involved in the eastward extrusion of the Tibetan Plateau [27], is a seismically active block [28]. Many strong earthquakes (M W > 6.5) in China occurred in this block and its border faults in past twenty years (Figure 1), including the 2001 M W 7.9 Kokoxili earthquake, the 2008 M W 7.8 Wenchuan earthquake, the 2010 Yushu M W 6.9 earthquake, the 2013 M W 6.6 Lushan earthquake, and the 2017 M W 6.5 Jiuzhaigou earthquake. In order to investigate geomagnetic anomalies before, during, and after earthquakes, the China Earthquake Administration operates a geomagnetic network including geomagnetic observatories and geomagnetic repeat stations in the Tibetan Plateau. However, only a repeat station located 26 km from the epicenter of the M W 6.6 Lushan earthquake recorded PMs [29]. PMs can be neglected if the fault rupture is small [30], unless observations are conducted right on faults [3]. However, PMs must be investigated carefully for earthquakes with ruptures larger than 100 km for their huge energy release and strong impact on local and regional geomagnetic fields.
in this block and its border faults in past twenty years (Figure 1), including the 2001 MW 7.9 Kokoxili earthquake, the 2008 MW 7.8 Wenchuan earthquake, the 2010 Yushu MW 6.9 earthquake, the 2013 MW 6.6 Lushan earthquake, and the 2017 MW 6.5 Jiuzhaigou earthquake. In order to investigate geomagnetic anomalies before, during, and after earthquakes, the China Earthquake Administration operates a geomagnetic network including geomagnetic observatories and geomagnetic repeat stations in the Tibetan Plateau. However, only a repeat station located 26 km from the epicenter of the MW 6.6 Lushan earthquake recorded PMs [29]. PMs can be neglected if the fault rupture is small [30], unless observations are conducted right on faults [3]. However, PMs must be investigated carefully for earthquakes with ruptures larger than 100 km for their huge energy release and strong impact on local and regional geomagnetic fields. The MW 7. 3 Maduo earthquake (herein referred to as Maduo earthquake) occurred at 18:04 UTC on 21 May 2021, with the epicenter located at 98.34° E and 34.59° N (http://www.csi.ac.cn/, accessed on 10 November, 2021). Maduo earthquake caused an about 160 km rupture ( Figure 1) along the Maduo-Gande fault [31]; however, the tectonic deformation is significantly weak in the region [32]. How does the earthquake disturb geomagnetic fields in the Bayan Har block? Are the piezomagnetic fields localized? More importantly, how will the earthquake impact observations of piezomagnetic fields in the future? To address these questions, in this study, we construct a piezomagnetic model to calculate the static coseismic PMs based on a slip model of the Maduo earthquake. Moreover, observed data of 2019 and 2021 on a repeat station are used to examine the PMs.   [31]; however, the tectonic deformation is significantly weak in the region [32]. How does the earthquake disturb geomagnetic fields in the Bayan Har block? Are the piezomagnetic fields localized? More importantly, how will the earthquake impact observations of piezomagnetic fields in the future? To address these questions, in this study, we construct a piezomagnetic model to calculate the static coseismic PMs based on a slip model of the Maduo earthquake. Moreover, observed data of 2019 and 2021 on a repeat station are used to examine the PMs.

Data and Methods
There are two strategies to calculate the PMs due to earthquake rupture. One is Stacey's scheme, involving stress-magnetic calculation [2]. Another is Green's function method, involving strain-magnetic calculation [33]. In this study, we prefer Stacey's scheme, because the coseismic stress change is easily calculated nowadays. A piezomagnetic model is constructed based on a slip model of Maduo earthquake to calculate the PMs.

Data and Methods
There are two strategies to calculate the PMs due to earthquake rupture. One is Stacey's scheme, involving stress-magnetic calculation [2]. Another is Green's function method, involving strain-magnetic calculation [33]. In this study, we prefer Stacey's scheme, because the coseismic stress change is easily calculated nowadays. A piezomagnetic model is constructed based on a slip model of Maduo earthquake to calculate the PMs.

Coseismic Slip Model
Zhao et al. determined the slip distribution of the coseismic rupture of Maduo earthquake based on InSAR data ( Figure 2) [32]. Their results show that the length of surface rupture is about 160 km distributed on five sub-vertical segments with strikes of 89.5°, 106°, 97.6°, 117.5°, and 83° and dips of 90°, 75°, 80°, 90°, and 90°, respectively. The faulting style is primarily left-lateral. The maximum strike slip and total slip (the sum of strike slip and dip slip) are about 6 m and 8 m at 6 km depth on segment 4, respectively. Most of them are released in the upper crust (0-15 km). There are three areas with significant coseismic slip (larger than 5 m), two areas in segment 2, and the other one is near the triple junction zone between segment 3, segment 4, and segment 5 ( Figure 2).

Layer Elastic Model
We adopt the PSGRN/PSCMP code to calculate the coseismic stress changes based on an elastic crust model [34]. The crust model is based on the results of a previous study [35]. The detailed parameters of the model are listed in Table 1. The crust is discretized into 3864 sub-blocks with a dimension of 2 × 2 × 2 km.

Layer Elastic Model
We adopt the PSGRN/PSCMP code to calculate the coseismic stress changes based on an elastic crust model [34]. The crust model is based on the results of a previous study [35]. The detailed parameters of the model are listed in Table 1. The crust is discretized into 3864 sub-blocks with a dimension of 2 × 2 × 2 km.

Piezomagnetic Model
The constitutive law of the stress and magnetization has been proposed as follows [36]: where ∆J is the change of magnetization, S is the deviatoric stress tensor, J is the initial magnetization, and β is the stress sensitivity. The piezomagnetic potential dW at point P (x 0 , y 0 , z 0 ) due to ∆J of a sub-block with coordinates of (x, y, z) can be expressed as Equation (2) [36], and the PMs can be obtained by differentiating the right-hand side of Equation (2).
where (l, m, n) are the direction cosines of J, and ρ is the distance between point P and the sub-block, The total PMs at point P are equal to the sum of the PMs derived from each sub-block, as follows: where H is calculation depth, The total intensity of PMs can be obtained through declination (D) and inclination (I) as follows [36]: The direction of initial magnetization is assumed to be consistent with the present magnetic fields, which have an inclination and a declination of 54 • and −1 • (westward), respectively. The magnetic anomalies are not strong in the Tibetan Plateau [37], which indicates a weak heterogeneity in magnetization. Thus, a uniform initial magnetization seems reasonable. Stress sensitivity of the order of 10 −3 MPa −1 is considered in many studies of seismomagnetic anomalies [3,7,13,19]. However, in some studies, the initial magnetization of 1A/m and stress sensitivity of 1 × 10 −3 MPa −1 are too small to explain the observed data [21,25,29]. In our model, an average magnetization of 3 A/m and a stress sensitivity of 2 × 10 −3 MPa −1 are assumed. A depth of 16 km was chosen because most of the stress associated with Maduo earthquake released in the upper crust (0-15 km).

Results
In order to eliminate the distortion of PMs resulting from stress concentration and singularity on the fault plane, we do not consider the PMs in the project area of fault slip about 2 km from the fault trace. The PMs are calculated at 1 m, where magnetometers are usually installed. The total intensity (F) and east component (Y) of PMs are calculated because Y is expected to be significantly changed under the east-west slip regime.
The spatial distribution of the PMs is shown in Figure 3. Our results show that the pattern of PMs is consistent with the left-lateral slip. For Y, negative changes appear in the compressional area, whereas positive changes appear in the extensional area. This is the opposite for F, because F is primarily calculated from the north component (X) and the vertical component (Z) which are reverse variations with Y according to Equation (1).

Results
In order to eliminate the distortion of PMs resulting from stress concentration and singularity on the fault plane, we do not consider the PMs in the project area of fault slip about 2 km from the fault trace. The PMs are calculated at 1 m, where magnetometers are usually installed. The total intensity (F) and east component (Y) of PMs are calculated because Y is expected to be significantly changed under the east-west slip regime.
The spatial distribution of the PMs is shown in Figure 3. Our results show that the pattern of PMs is consistent with the left-lateral slip. For Y, negative changes appear in the compressional area, whereas positive changes appear in the extensional area. This is the opposite for F, because F is primarily calculated from the north component (X) and the vertical component (Z) which are reverse variations with Y according to Equation (1).
The large PMs occur in the near-field. The maximum amplitude of PMs is about 10 nanotelsa (nT) for F and 8 nT for Y. Additionally, three peaks are in the regions with much larger slip mentioned above. However, PMs are small at the epicenter area. The large PMs occur in the near-field. The maximum amplitude of PMs is about 10 nanotelsa (nT) for F and 8 nT for Y. Additionally, three peaks are in the regions with much larger slip mentioned above. However, PMs are small at the epicenter area.
We also examine the Y changes on three profiles (Figure 4). Since the distributions of the F changes are opposite with the Y changes, we do not further analyze the F changes. The Y changes are more remarkable at the north side than that at the south side of fault trace. The Y changes on profile 1 and profile 2 are about zero at 20 km and 30 km from the fault trace on the south side and north side, respectively. The Y changes on the profile 3 are about 0.5 nT at 40 km away from the fault trace.
We also examine the Y changes on three profiles (Figure 4). Since the distributions of the F changes are opposite with the Y changes, we do not further analyze the F changes. The Y changes are more remarkable at the north side than that at the south side of fault trace. The Y changes on profile 1 and profile 2 are about zero at 20 km and 30 km from the fault trace on the south side and north side, respectively. The Y changes on the profile 3 are about 0.5 nT at 40 km away from the fault trace.

Validation of Piezomagnetic Model
We compare the PMs of the Maduo earthquake with that of other major continental strike-slip earthquakes (M > 6), including the 1989 ML 7.1 Loma Prieta earthquake [19], the 1992 MW 7.3 Landers earthquake [20], and the 2004 M 6.0 Parkfield earthquake [4]. It should be noticed that the PMs strongly depend on parameters. Therefore, we focus on relative values of PMs here.
The PMs are proportional to initial magnetization and stress sensitivity according to Equation (1). We calculate the normalized PMs (NPMs) as follows: The PMs exhibit a negative correlation with depth of slip (Hs) [38] and positive correlation with slip (S). The earthquake magnitude or moment considered by Mueller and Johnston [29] is substituted by the ratio of S and Hs in this paper. Because we cannot obtain the value of stress sensitivity for the Loma Prieta earthquake, 1 × 10 −3 MPa −1 and 2 × 10 −3 MPa −1 are assumed in the piezomagnetic model of Loma Prieta. Here, we choose the maximum PMs of these earthquakes to analyze. The correlation between maximum NPMs and S/Hs is shown in Figure 5. All values of these piezomagnetic models are shown in Table 2. Figure 5 shows that NPMs increase with S/Hs. In other words, given the same β and J, the larger S/Hs is, the larger the PMs are. As shown in Figure 5, NPMs of Maduo earthquake are consistent with those of other strike-slip earthquakes and the Loma Prieta earthquake with stress sensitivity of 2 × 10 −3 MPa −1 . It indicates that our piezomagnetic model is reasonable.

Validation of Piezomagnetic Model
We compare the PMs of the Maduo earthquake with that of other major continental strike-slip earthquakes (M > 6), including the 1989 M L 7.1 Loma Prieta earthquake [19], the 1992 M W 7.3 Landers earthquake [20], and the 2004 M 6.0 Parkfield earthquake [4]. It should be noticed that the PMs strongly depend on parameters. Therefore, we focus on relative values of PMs here.
The PMs are proportional to initial magnetization and stress sensitivity according to Equation (1). We calculate the normalized PMs (NPMs) as follows: The PMs exhibit a negative correlation with depth of slip (Hs) [38] and positive correlation with slip (S). The earthquake magnitude or moment considered by Mueller and Johnston [29] is substituted by the ratio of S and Hs in this paper. Because we cannot obtain the value of stress sensitivity for the Loma Prieta earthquake, 1 × 10 −3 MPa −1 and 2 × 10 −3 MPa −1 are assumed in the piezomagnetic model of Loma Prieta. Here, we choose the maximum PMs of these earthquakes to analyze. The correlation between maximum NPMs and S/Hs is shown in Figure 5. All values of these piezomagnetic models are shown in Table 2. Figure 5 shows that NPMs increase with S/Hs. In other words, given the same β and J, the larger S/Hs is, the larger the PMs are. As shown in Figure 5, NPMs of Maduo earthquake are consistent with those of other strike-slip earthquakes and the Loma Prieta earthquake with stress sensitivity of 2 × 10 −3 MPa −1 . It indicates that our piezomagnetic model is reasonable.

Possible Values of Initial Magnetization and Stress Sensitivity
Undoubtedly, we cannot obtain a completely correct estimation of PMs using uniform parameters in a piezomagnetic model [12]. The precise structure of magnetization is not easy to determine. Nevertheless, we still attempt to evaluate the initial magnetization and stress sensitivity assumed in our piezomagnetic model.
We examined the observed geomagnetic data on a geomagnetic repeat station (MD) located about 6 km from the fault trace, as shown in Figure 1. A non-magnetic plastic plate was set on the repeat station to ensure that the magnetometers were mounted on the same position every time. The preseismic observation was conducted at 06:09 to 06:41 UTC on 19 June 2019 to investigate the regional geomagnetic fields. The repeat observation was conducted at 06:05 to 06:31 UTC on 27 June 2021 to investigate the geomagnetic field changes associated with the Maduo earthquake. The Kp indices (the planetary 3-h-range Kp index was designed to measure solar particle radiation by its magnetic effects) were both 1 at two observations times, respectively, which suggests that there were no strong magnetic disturbances, solar flares, or solar storms. The geomagnetic vector data were obtained including declination (D), inclination (I), and total intensity (F). GSM-19T proton magnetometer with an accuracy of 0.2 nT was used to measure F. Flux-gate theodolite with an accuracy of 0.2 arcmin was used to measure D and I. Y was calculated according to D, I, and F, as Y = F cos I sin D. In order to calculate the geomagnetic field changes, we must eliminate the secular and diurnal variations, which are more than a few tens of nanotesla [39]. The methods introduced by Wang et al. [40] including diurnal reduction and secular reduction are used to correct the raw data. The changes of F and Y between 2019 and 2021 are −0.55 nT and +4.8 nT, respectively, which seem to be large enough to be detected by magnetometers.
Here, we must investigate the error of repeat observation by analyzing the mounting error of magnetometers, diurnal reduction error, and secular reduction error. The maximum magnetic gradient of MD is 0.2 nT/m. Thus, the mounting error resulting from position deformation can be neglected. The larger RMS of diurnal reduction in 2019 and 2021

Possible Values of Initial Magnetization and Stress Sensitivity
Undoubtedly, we cannot obtain a completely correct estimation of PMs using uniform parameters in a piezomagnetic model [12]. The precise structure of magnetization is not easy to determine. Nevertheless, we still attempt to evaluate the initial magnetization and stress sensitivity assumed in our piezomagnetic model.
We examined the observed geomagnetic data on a geomagnetic repeat station (MD) located about 6 km from the fault trace, as shown in Figure 1. A non-magnetic plastic plate was set on the repeat station to ensure that the magnetometers were mounted on the same position every time. The preseismic observation was conducted at 06:09 to 06:41 UTC on 19 June 2019 to investigate the regional geomagnetic fields. The repeat observation was conducted at 06:05 to 06:31 UTC on 27 June 2021 to investigate the geomagnetic field changes associated with the Maduo earthquake. The Kp indices (the planetary 3-h-range Kp index was designed to measure solar particle radiation by its magnetic effects) were both 1 at two observations times, respectively, which suggests that there were no strong magnetic disturbances, solar flares, or solar storms. The geomagnetic vector data were obtained including declination (D), inclination (I), and total intensity (F). GSM-19T proton magnetometer with an accuracy of 0.2 nT was used to measure F. Flux-gate theodolite with an accuracy of 0.2 arcmin was used to measure D and I. Y was calculated according to D, I, and F, as Y = F cos I sin D. In order to calculate the geomagnetic field changes, we must eliminate the secular and diurnal variations, which are more than a few tens of nanotesla [39]. The methods introduced by Wang et al. [40] including diurnal reduction and secular reduction are used to correct the raw data. The changes of F and Y between 2019 and 2021 are −0.55 nT and +4.8 nT, respectively, which seem to be large enough to be detected by magnetometers.
Here, we must investigate the error of repeat observation by analyzing the mounting error of magnetometers, diurnal reduction error, and secular reduction error. The maximum magnetic gradient of MD is 0.2 nT/m. Thus, the mounting error resulting from position deformation can be neglected. The larger RMS of diurnal reduction in 2019 and 2021 observed data are 0.11 arcmin, 0.03 arcmin, and 0.18 nT for D, I, and F, respectively, which can be regarded as the error of diurnal reduction. The sixth-order natural orthogonal component (NOC) of the geomagnetic fields from June 2019 to July 2021 is used to correct the secular variation. Figure 6 depicts the observed data in 72 quietest days and the values of NOC model from June 2019 to July 2021 on Jiayuguan geomagnetic observatory, which is one of the observatories used in the NOC model. The RMS of difference between observed data and NOC model is 0.08 arcmin, 0.1 arcmin, and 0.75 nT for D, I, and F, respectively, which can be regarded as the error of secular reduction. The total error can be calculated by the accuracy of instruments, diurnal reduction error, and secular reduction error, which 8 of 11 are 0.24 arcmin, 0.22 arcmin, and 0.80 nT for D, I, and F, respectively. As Y is primarily dependent on D, the total error of Y can represent by the error of D, which is about 2.2 nT in MD for 0.24 arcmin.
which are 0.24 arcmin, 0.22 arcmin, and 0.80 nT for D, I, and F, respectively. As Y is primarily dependent on D, the total error of Y can represent by the error of D, which is about 2.2 nT in MD for 0.24 arcmin.
Since the observed F (−0.55 nT) is less than the error (0.80 nT), we do not compare it with calculated F changes. As the observed Y is about twice the error, we believe it was reasonable and authentic. The observed and predicted Y at MD are +4.8 ± 2.2 nT and +4.3 nT, respectively. We consider that our piezomagnetic model can explain the observed data. In other words, the observed geomagnetic changes before and after Maduo earthquake at MD are caused by stress changes resulting from fault rupture. It also indicates that initial magnetization of 3 A/m and stress intensity of 2 × 10 −3 MPa −1 assumed in piezomagnetic model of Maduo earthquake are reasonable.  Since the observed F (−0.55 nT) is less than the error (0.80 nT), we do not compare it with calculated F changes. As the observed Y is about twice the error, we believe it was reasonable and authentic. The observed and predicted Y at MD are +4.8 ± 2.2 nT and +4.3 nT, respectively. We consider that our piezomagnetic model can explain the observed data. In other words, the observed geomagnetic changes before and after Maduo earthquake at MD are caused by stress changes resulting from fault rupture. It also indicates that initial magnetization of 3 A/m and stress intensity of 2 × 10 −3 MPa −1 assumed in piezomagnetic model of Maduo earthquake are reasonable.

Potential Usefulness of Observations along Fault
Many magnetometers fail to record PMs during earthquakes. On one hand, they are far away from fault rupture. As our piezomagnetic model shows, even for a M W 7.3 earthquake with a rupture length of 160 km, the geomagnetic field changes that are above the magnetometer accuracy are limited in areas 40 or 50 km from rupture. Considering the error resulting from field observations and data processing, we consider that only magnetometers mounted at a region within 20 km from the fault rupture can observe useful seismomagnetic signals. On the other hand, it also could fail to record seismomagnetic signals even when the magnetometers are mounted in near-field. As our results show, the PMs are small in epicenter area. Thus, only stations that are ideally located can observe the PMs.
The near-field PMs show variations in the degree of localized strengthening along rupture trace (Figure 3), especially at the triple junction of three fault segments. Indeed, the PMs are expected to be larger at ends of faults or locations with geometrical complexities, such as fault bend or branching [40]. Bifurcation of slip appears to be a common feature of some major continental strike-slip earthquakes [41,42]. It indicates that locations with geometrical complexities are potential sites to observe the PMs.

Conclusions
The 2021 M W 7.3 Maduo earthquake is assumed to disturb local and regional geomagnetic fields significantly. A piezomagnetic model is constructed to investigate the geomagnetic fields changes due to stress changes resulted from earthquake rupture. We validate the piezomagnetic model by comparing with piezomagnetic changes of other strike-slip earthquakes and observed data before and after earthquake on a geomagnetic repeat station. The PMs of Maduo earthquake seem to appear in areas within 40 km from fault trace. The PMs are larger at the location with bifurcation of slip. The observed geomagnetic field changes before and after Maduo earthquake can be explained by our piezomagnetic model with initial magnetization of 3 A/m and stress sensitivity of 2 × 10 −3 MPa −1 . Thus, if there are more observations, we can constrain the detailed magnetic properties and the slip model.