Determination of the local magnitudes of small earthquakes using a dense seismic array in the Changning–Zhaotong Shale Gas Field, Southern Sichuan Basin

With the development of unconventional shale gas in the southern Sichuan Basin, seismicity in the region has increased significantly in recent years. Though the existing sparse regional seismic stations can capture most earthquakes with ML≥2.5 , a great number of smaller earthquakes are often omitted due to limited detection capacity. With the advent of portable seismic nodes, many dense arrays for monitoring seismicity in the unconventional oil and gas fields have been deployed, and the magnitudes of those earthquakes are key to understand the local fault reactivation and seismic potentials. However, the current national standard for determining the local magnitudes was not specifically designed for monitoring stations in close proximity, utilizing a calibration function with a minimal resolution of 5 km in the epicentral distance. That is, the current national standard tends to overestimate the local magnitudes for stations within short epicentral distances, and can result in discrepancies for dense arrays. In this study, we propose a new local magnitude formula which corrects the overestimated magnitudes for shorter distances, yielding accurate event magnitudes for small earthquakes in the Changning−Zhaotong shale gas field in the southern Sichuan Basin, monitored by dense seismic arrays in close proximity. The formula is used to determine the local magnitudes of 7,500 events monitored by a two‐phased dense array with several hundred 5 Hz 3C nodes deployed from the end of February 2019 to early May 2019 in the Changning−Zhaotong shale gas field. The magnitude of completeness ( MC ) using the dense array is −0.1, compared to MC 1.1 by the sparser Chinese Seismic Network (CSN). In addition, using a machine learning detection and picking procedure, we successfully identify and process some 14,000 earthquakes from the continuous waveforms, a ten‐fold increase over the catalog recorded by CSN for the same period, and the MC is further reduced to −0.3 from −0.1 compared to the catalog obtained via manual processing using the same dense array. The proposed local magnitude formula can be adopted for calculating accurate local magnitudes of future earthquakes using dense arrays in the shale gas fields of the Sichuan Basin. This will help to better characterize the local seismic risks and potentials.


Introduction
Earthquake magnitude is a fundamental parameter for characterizing earthquake sources and estimating source energy, seismic moment, source rupture time, rupture size and area (Yan ZD, 1992). The original magnitude scale proposed by Richter (1935) for Southern California is commonly referred as the local magnitude , which is widely used in various countries. Li SB (1981) introduced the local magnitude into China in 1959 and summarized the calibration functions suitable for seismometers of Type 62 (short period seismograph) and Type SK (medium and long period seismograph) commonly used in China during that time. The function describes the attenuation characteristics of seismic waves with epicentral distance and is closely related to the crustal structure (Chen PS and Qin JZ, 1983;Xue ZZ, 1992). It should be noted that the calibration function does not consider the influence of the focal depth, since most of the earthquakes in China are shallow intraplate continental earthquakes and the hypocentral depth is usually negligible compared to the epicentral distance (Liu RF et al., 2015). With the expansion of the domestic Chinese Seismic Network (CSN) (Dai GH and An YR, 2020), Wang LY et al. (2016) profiled the seismic phase data from 1,308 stations dated from 1973 to 2002 in 31 provincial centers, recalculating the calibration function which was then adopted in the Chinese national standard (GB 17740−2017). However, the current national standard for determining the local magnitude was initially proposed for permanent seismic networks at regional scales. Since the calibration function uses a minimum step of 5 km, the resolution could be poor for characterizing the local magnitudes of events monitored by seismic arrays at distances ranging from few hundred meters up to 30 km (e.g., Meng LY et al., 2019).  (Zou CN et al., 2015). Especially with the accelerated shale gas development in the Sichuan Basin since 2011, seismic activity has increased dramatically, including a series of destructive earthquakes in Junlian, Xingwen, Gongxian, and Changning (Yi GX et al., 2019, Lei XL et al., 2020 counties. Among them, the Xingwen earthquake occurring on December 16, 2018 has set a new record for induced earthquakes related to hydraulic fracturing (Atkinson et al., 2020). In contrast, there were much fewer felt earthquakes in the southern Sichuan Basin in the last decade while small earthquakes were continuously monitored by the national networks operated by the China Earthquake Administration (CEA) (Dai GH and An YR, 2020), with the magnitude of completeness ( ) being (Lei XL et al., 2019a;Meng LY et al., 2019). To better address the concern for increasing seismic activity, Meng LY et al. (2019) deployed six temporary seismometers with about 20 km interstation spacing covering the Gongxian county in the Changning-Zhaotong shale gas field between February 2015 and April 2017. In May 2017, with the deployment of 15 additional seismometers in the same area, a regional seismic network consisting of 21 stations started to monitor the Changning area, including the Changning shale gas field and the adjacent Changing salt mine. The average interstation spacing of the expanded regional network is about 12 km, enabling detection of smaller earthquakes. From 2015 to 2017, more than 15,000 earthquakes with were monitored by this regional network, about 86% of which were located in the Changning−Zhaotong shale gas field. The of those events is 1.1, much smaller than those reported in the China National Network Catalog (CNNC) (Dai GH and An YR, 2020).
With the advent of portable seismic nodes, dense seismic arrays have been widely used to evaluate the seismic risks and potentials of unconventional shale gas and oil development through M L microearthquake monitoring (e.g. Bommer et al., 2006;Majer et al., 2012). To monitor the seismicity likely related to hydraulic fracturing (Lei XL et al., 2019b and characterize reactivation of the fault system in the Changning-Zhaotong shale gas field, we deployed a two-phased Dense Seismic Array (DSA) with about 1.5 km interstation spacing on average between 28th February 2019 and 6th May 2019, lasting about 70 days. The epicentral distances of this dense array range from several hundred meters to about 20 km. The amplitude of the body-wave decays rapidly with epicentral distance in the first few kilometers due to a combination of contributions from geometrical spreading, attenuation, source radiation pattern, etc. Thus, if the traditional local magnitude formula is used, the individually determined local magnitudes using the body-wave amplitude at different stations in close proximity will vary rapidly, since the calibration function (with a minimal step of 5 km) is unable to compensate for the amplitude variation. Several studies also found overestimation in magnitudes at close distances in other surveys (Scognamiglio et al., 2012;Butcher et al., 2017;Luckett et al., 2019). For instance, Butcher et al. (2017) showed that this discrepancy is significant because the travel paths are mostly within the sedimentary layer rather than the underlying basement for a station nearby, and thus a different attenuation or calibration term is required in estimating the scale.
To accurately calculate the local magnitude for events in the Changning−Zhaotong shale gas field, we propose a new formula which is applicable for stations at varying epicentral distances. Many previous studies were primarily focused on verifying the southern California scale originally proposed by Richter (1935), or to recalibrate it for different areas to take into account diverse attenuation properties (Hutton and Boore, 1987;Langston et al., 1998;Keir et al., 2006;Di Bona, 2016). In our study, hundreds of common events observed by both the DSA and CSN are used to derive the geometric spreading and attenuation coefficients in the calculation formula using a least-squares method, where the magnitudes provided by the CNNC are regarded as the reference. The CNNC is selected as the reference because 1) the magnitudes for many small events in the Changning−Zhaotong shale gas field are only available in this catalog, 2) the magnitudes in CNNC are considered as the authoritative result given by the China Earthquake Networks Center, 3) the seismic records are more continuous, and 4) the calculation method for the magnitudes is uniform.
In addition, the rapid growth in seismic data volume poses increasing burden on processors, and many small events may be inevitably omitted (Wang RJ et al., 2020). To fully extract and characterize earthquakes from continuous waveform data and better delineate fault reactivations, we also adopt a machine-learning strategy for automated event detection and phase picking, followed by accurate event location and magnitude determination. The machine-learning based strategy not only shows comparable accuracy compared to the manual processing, but also detects nearly twice the events. Both the minimal and the of the detected events are further reduced.

Data
Between 28th February 2019 and 6th May 2019, we deployed a two-phased seismic dense array with ~1.5 km interstation spacing in the Changning−Zhaotong shale gas field in the southern Sichuan Basin, to investigate the frequent seismicity likely related to hydraulic fracturing (Lei XL et al., 2017, 2019b. The first phase of the dense array deployed from 28th February to 3rd April 2019 consists of 187 short-period seismic nodal stations, and the following second phase of the dense array deployed from 4th April to 6th May consists of 150 stations, covering an area of nearly 800 in total (Figure 1). Two types of three-component 5 Hz nodal stations, SmartSolo and Zland, were deployed in a mixed distribution. These two types of seismic nodes have very similar instrument responses, and the sampling frequency for all the nodes is set to 500 Hz or 2 ms. Since there is no available information for the local noise level, the expected microseismic magnitude, or other important factors for the passive seismic survey, the gain of the Zland nodes was set to , while that of the SmartSolo was set to to acquire data in a wider amplitude range and to accommodate various field situations (Fan TJ et al., 2019). Since the gain of the SmartSolo nodes was set at the maximum level, the waveform amplitudes at those nodes can be saturated for larger earthquakes. Thus, data from those nodes are excluded and only data from the Zland nodes with 12 dB gain are included when earthquakes with magnitude greater than are used to derive the local magnitude formula. j To detect the events, we cut the continuous waveform data from the dense array into one-minute segments, and in each segment the kurtosis value is evaluated for the Z-component of station (Saragiotis et al., 2002): x E μ σ where is a one-minute waveform, is the expectation, is the mean value and is the deviation. For a normal distribution (Gaussian noise), the Kurtosis value is 3, and for a trace containing an earthquake signal the value is larger than 3. However, spiky noise and other interference deviating from the normal process can also increase the kurtosis value of a trace and result in false detection. Thus, to mitigate false detection, the median of the kurtosis values of all traces in a one-minute segment is used, with a threshold set to 10 for our event detections. All the 9,758 initiallyclaimed detections by the kurtosis value method are further checked manually, and only a few hundred are false claims. Eventually, 7,543 events are manually processed, yielding 506,652 P-arrival phases and 433,639 S-arrival phases. We use DSA-M (Dense Seismic Array-Manual) to refer to this catalog in the following discussions. In comparison, 1,292 seismic events ( Figure 2) occurring in the study area have been cataloged in the CNNC for the same 70-day period. Figure 2 shows the locations of the earthquakes detected by our DSA and provided by the CNNC using the double-difference method (Zhang HJ andThurber, 2003, 2006) based on a local velocity model obtained from surface wave ambient noise tomography using data from the same array (Li JL et al., 2019). The inset shows the P-and S-wave traveltime residuals, both of which are in Gaussian distributions. The standard devi-  , and the yellow one shows the dense network of the second phase (phase-2). The black triangles denote the Zland nodes and the blue triangles indicate the SmartSolo nodes. The cyan triangles denote the regional stations deployed by Meng LY et al. (2019) since 2015. The Changning (north) and the Zhaotong (south) shale gas development zones are indicated by the irregularly shape areas marked by dashed white lines. The dense array covered two previous major earthquakes, the Xingwen earthquake (2018/12/16, ) and the Gongxian earthquake (2019/1/3 ).

National Standard for Local Magnitude (M L ) Calculation in China
The magnitude scale was first defined by Richter (1935) for Southern California, and all subsequent magnitude scales are related to this one (Havskov and Ottemoller, 2010). Since the CSN did not install the Wood-Anderson standard seismometer in the 1950s, the scale proposed by Richter could not be directly transferred to China. In 1959, the local magnitude was introduced into China (Li SB, 1981) based on seismometer Type 62 (short period seismograph) and seismometer SK (medium and long period seismograph), which were commonly used in China in the 1950s (Liu RF et al., 2015). The local magnitude measured at a station is defined as: and are the calibration function and the epicentral distance in kilometers, respectively, and is the arithmetic mean of the maximum S-wave amplitudes at the two horizontal components in micrometers: where is the maximum S-wave or Lg-wave amplitude in the north-south component in μm and is the maximum S-wave or Lg-wave amplitude in the east-west component in μm. The final local magnitude for an event is the arithmetic mean of all the magnitudes measured at individual stations. According to the Chinese national standard (GB 17740−2017) for earthquake magnitude, the seismic record should be simulated to the short period displacement record of a DD-1 seismometer before the magnitude calculation, and the maximum amplitude should be more than twice the interference level. The transfer function of the DD-1 short-period seismograph is: , and is the frequency.
The calibration function in Equation (2) describes the attenuation characteristics of seismic waves with epicentral distance, which is closely related to the propagation path (Chen PS and Qin JZ, 1983;Xue ZZ, 1992). For stations deployed in close proximity to the event where the seismic amplitude decays rapidly, a refined local magnitude calibration function should be established, since an inaccurate calibration function may lead to varying measured magnitudes for the same earthquake at different stations (Chen PS and Qin JZ, 1983). Recently, Wang LY et al. (2016) derived the regional calibration function which was then adopted in the Chinese national standard (GB 17740−2017) with seismic phase data from 1,308 stations in 31 provincial centers dating from 1937 to 2002, including 105,282 effective earthquakes above based on the statistical analysis for magnitude residuals (Chen PS and Qin JZ, 1983). The 31 provinces in China are divided into five regions, each of which has its own adjusted calibration function. Those regions are: Northeast-North China (Heilongjiang, Jilin, Liaoning, Inner Mongolia, Beijing, Tianjin, Hebei, Shanxi, Shandong, Ningxia, and Shaanxi), South China (Fujian, Guangdong, Guangxi, Hainan, Jiangsu, Shanghai, Zhejiang, Jiangxi, Hunan, Hubei, and Anhui), Southwest China (Yunnan, Sichuan, Chongqing, and Guizhou), Qinghai−Xizang (Qinghai, Tibet, and Gansu) and Xinjiang. In the updated calibration function, the epicentral distance ranges from 0 to 1000 km. Since the maximum epicentral distance for our two-phased dense array in the Changning−Zhaotong shale gas field is less than 40 km, only the relevant calibration function is listed in the following Table 1. Figure 3 shows 1,223,274 epicenter distances from the 7,543 earthquakes recorded by the two-phased dense array. Figure 3a shows that if the same 5 km minimal step as in the calibration function from the national standard (GB 17740−2017) is used (Table 1), the varying epicentral distances can only be divided into seven different groups, while most epicentral distances fall into four of those. In comparison, when the epicentral distances are further divided with 1 km or 0.5 km refined intervals, it is found the variation in epicentral distances is much more complicated. The comparison indicates that the currently adopted calibration function with large intervals cannot accurately reflect the variation in epicentral distance for stations with relatively short epicentral distances, and can result in biased magnitude calculations at stations in close proximity since the body-wave amplitude decays rapidly.
We further elaborate on this issue by calculating the local magnitude of an earthquake monitored by the dense array using Equation (2). To calculate the local magnitude at a given station, we first remove the instrument response of the short period nodes (available at the IRIS website, http://ds.iris.edu/NRL/), and then convert the velocity record into displacement, followed by convolving the displacement with the transfer function of DD-1 in Equation (4). The maximum amplitude in a 3 second window around the S-wave arrival time is then used for Equation (3). The magnitudes at various stations where we manually pick the Swave arrivals (colored triangles) are shown in Figure 4. Three circles with radii of 5 km, 10 km, 15 km, respectively, are shown with gray dashed lines. However, since the values of the calibration function are 2.0 for all the distances (Table 1), the calculated magnitudes at different stations in the dense array vary considerably. The reported magnitude in the CNNC using more distant stations for this particular event is , which is close to the magnitude estimated by our stations situated at the outer circles around .
The lower right inset also shows the magnitude differences against the epicentral distances, which clearly indicates that the stations within 10 km tend to overestimate the magnitude in general. The estimated by averaging the local magnitudes measured at various stations (colored triangles) is 2.4, which is a significant overestimate.
In addition, we also analyze the magnitude discrepancy ( ) against the epicentral distance at different stations of the dense array for the 995 common events in the DSA-M and the CNNC catalogs ( Figure 5). Each square indicates the average of the magnitude discrepancy for all the events at the same epicentral distances, and the error bar shows one standard deviation. The local magnitudes at different stations in Figure 5 are calculated by Equation (2). It is clear that there is a significant overestimate in magnitude at stations less than 10 km (up to oneunit of magnitude), whereas there is a systematic underestimate in magnitude at stations further than 12 km. Figures 3, 4 and 5 indicate that the calibration function with a minimal step of 5 km is not suitable for characterizing the local magnitude of events monitored by dense seismic array, whose epicentral distances usually vary from a few hundred meters to less than 30 km.

Magnitude Calculation with the Dense Seismic Array (DSA)
Since most of the earthquakes in China are shallow intraplate continental earthquakes (Liu RF et al., 2015), the current calibration      the magnitude measurement. Thus, the hypocentral distance is used instead in the following calculations.
The original formula for calculating the local magnitude is defined as (e.g., Havskov and Ottemoller, 2010): where is the zero-to-peak displacement amplitude at a given station with epicentral distance . Here, we use the simulated Swave displacement amplitude as measured by the DD-1 seismometer.
is the displacement correction term, i.e., the calibration function. The amplitude of the S-wave can be expressed as a function of hypocentral distance : where is the initial amplitude at the source, is the geometric spreading factor, which is 1 for body waves and 0.5 for surface waves, is the frequency, is the S-wave velocity, and is the quality factor. It may not be easy to directly determine the values of those factors, and thus empirical estimations are oftentimes used. Taking logarithm on both sides of Equation (6) yields Here, the distance-related amplitude attenuation term is defined as: (5), and if and are assumed to be constants, then can be written as: where are constants representing geometrical spreading, attenuation, and the base level, respectively. Note that the values of those parameters depend on the particular survey region and should be adjusted accordingly. In the following, we invert for proper coefficients in Equation (9) and propose a new magnitude calculation formula suitable for stations with short epicentral distances in the study area. As mentioned above, the local magnitude for event is determined by taking the arithmetic average of magnitudes measured at each individual station . That is, for the earthquake, the local magnitude measured at the station is: M i and local magnitude for this event is: or: where is the number of stations that participated in calculating the event magnitude .
In the study area, there are 995 relatively strong, common earthquakes recorded by both the CNNC and DSA-M between February 28th 2019 and May 6th 2020. The magnitudes of these earthquakes determined by our dense array with shorter hypocentral distances should not be altered with the new formula, i.e. magnitudes provided by the CNNC should be used as the reference. A least-squares fitting is used to derive the coefficients in Equation (12) by minimizing the magnitude difference between from our DSA and the reference from the CNNC. The objective function for the magnitude difference can then be expressed as: where is the number of common events in both the CNNC and DSA-M. Optimal solutions for the coefficients in Equation (13) can be found by minimizing the objective function : Rearranging Equation (14), then we get: The coefficients can be obtained by solving Equation (15).
For a study area covered by a dense seismic array in China, the procedure for estimating the local magnitude can be summarized in the following: (1) Find the common events in the CNNC catalog provided by China Earthquake Administration (CEA) and also recorded by the dense array; (2) Preprocess the seismic waveform (e.g. remove the mean and linear trend, taper the head and tail of a waveform, remove instrument response, convert velocity into displacement) and simulate recording to a DD-1 short-period seismometer by convolving with the DD-1 instrument response (Equation (4)); m 1 , m 2 , m 3 (3) Optimize the parameters by solving Equation (15); (4) With the best-determined coefficients , process the station waveform with step (2) and calculate the magnitude by Equation (12). Figure 6 shows the detailed flowchart for estimating the local magnitudes of earthquakes monitored by a dense array.
Using the 995 common earthquakes with magnitudes referenced to the CNNC, the local magnitude formula for the Changning− Zhaotong shale gas field in the southern Sichuan Basin is: To further validate the proposed new magnitude formula in Equation (16), we analyze the magnitude discrepancy (ΔM = M station -M CNNC ) against the hypocentral distances at different stations of the dense array for all the 995 common events in the DSA-M and the CNNC catalogs (Figure 7). Compared with Figure 5 which are calculated using Equation (2), it is obvious that the new formula successfully compensates for the variation of amplitude with distance, yielding relatively consistent magnitudes for stations at varying hypocentral distances. Appendix A compares the estimated moment and local magnitudes, and Appendix B further eval-uates the individual contributions from the near, intermediate and far-field terms of a point dislocation source.
Using the new magnitude formula, we estimate the relationship between earthquakes' occurrence and their magnitudes (Gutenberg-Richter law), which is: where N is the cumulative number of earthquakes with magnitude larger or equal to in a given period of time, defines the productivity, and is the earthquake frequency ratio of small to large events. In most tectonically active regions, the b-value is about (Scholz, 2015). In this study, we select the magnitude of completeness as the minimum magnitude when fitting the Gutenberg-Richter law to avoid bias and deviation in fitting caused by incomplete observation of the earthquakes. It is found that the distribution of the magnitudes of the 995 common earthquakes recorded by both DSA-M and the CNNC are rather similar, and the values of and in both catalogs are also quite close (Figure 8).
Moreover, we calculate the standard deviation for the local magnitudes in our DSA-M for the 995 common earthquakes according to Chen PS and Qin JZ (1983): The results are and using Equation (2) and Equation (16), respectively, indicating that it is more consistent and accurate to calculate earthquake magnitudes using the proposed formula for

Magnitude Calibration
Magnitude Calculation Figure 6. Flowchart summarizing the main steps in magnitude calibration and determination for DSA. The red dotted box shows the magnitude calibration procedure for DSA using a least-squares fitting. The black dotted box shows the procedure for calculating earthquake magnitudes. data acquired by dense seismic arrays deployed in close proximity.

Shale Gas Field
Event detection using the kurtosis values in general is effective for the data from the dense array. However, arrival-picking, which is rather time-consuming, is the prerequisite for magnitude calculation. Indeed, it took several months for us to manually pick the arrivals for thousands of earthquakes observed by the two-phased dense array, which yields 506,652 P-wave arrivals and 433,639 Swave arrivals in total. To further detect smaller events which are omitted in the detection using the kurtosis values, and reduce the subsequent efforts in arrival-picking, a machine-learning approach is used, which is briefly described herein.
First, Inception V3, a deep neural network for image recognition (Szegedy et al., 2015(Szegedy et al., , 2016 has been retrained for microseismic detection. Since three-component seismic data observed by a dense array resemble a RGB image, we use about 2,800 oneminute dense-array seismic data segments containing microseismic events and about 20,000 data segments without events to train a new top layer of Inception. The retrained Inception can be used to determine whether a given one-minute seismic data segment from the dense array contains one or more microseismic events or not. Using the retrained neural network specialized for microseismic event detection, we find about 14,000 small and me-dium earthquakes in total, nearly double the number of events detected by the kurtosis value method, and about a ten-fold increase compared to the CNNC in the study area over the same period. Subsequently, another neural network named PhaseNet (Zhu WQ and Beroza, 2019) is used to pick the P-and S-arrivals of the detected events. We refer to the catalog fully processed through machine-learning as DSA-AI in the following discussions. Admittedly, about 8.5% of the events in the DSA-M catalog which are detected by the kurtosis values and then manually processed are not included DSA-AI catalog. Most of those events are of small magnitudes with extremely low signal-to-noise ratio, and 78.7% of them have magnitudes less than 0.0. Table 2 lists the number of events in three different earthquake catalogs, DSA-M, DSA-AI and CNNC, and the number of common events among them. Finally, to verify the reliability of the automatic processing using machine learning, we compare the location results of the microearthquakes processed manually and automatically based on the same velocity model (Li JL et al., 2019) (Figure 9). The location results are highly consistent. Figure 10 shows the travel time-hypocenter distance curve of the 6,899 common events through the manual and   machine-learning based processing, which are also rather consistent.
Using Equation (16), we calculate the magnitudes of all the earthquakes in the DSA-M and DSA-AI ( Figure 11). The magnitude threshold in the catalog provided by CSN is , whereas the smallest earthquakes monitored by DSA processed manually and through machine learning are −1.0 and −1.3, respectively. Using the maximum curvature (MAXC) method (Mignan and Woessner, 2012), we find the magnitudes of completeness ( ) for the catalogs from the CNNC, DSA-M and DSA-AI are 1.1, −0.1 and −0.3, respectively, differing by more than one magnitude. Also, the b-value in the Gutenberg-Richter law (Bender, 1983) for the catalogs DSA-M and DSA-AI both are 1.01, practically the same with the b-value estimated from the CNNC.

Discussion
The number of events processed using machine learning (DSA-AI) almost doubles that processed manually (DSA-M). Small earth- quakes are key to better delineating the fault geometry (Fukuyama et al., 2003;Ross et al., 2017), the chain of foreshocks, mainshocks, and aftershocks (Shelly et al., 2016), and the triggering and nucleation of earthquakes (Cochran et al., 2004). The magnitudes for common events in DSA-M and DSA-AI, and for the extra events only in the DSA-AI are shown in Figure 12. When the magnitude is above 1.0, only a few events are omitted by the human processors compared to the machine-learning approach. However, when the event magnitude gets lower, an increasing amount of events are omitted by the human processors. Indeed, it seems to be difficult for human processors to further lower the processing threshold than −0.3, and much fewer small magnitude events can be processed in practice, since those events may be overlooked due to signal-to-noise ratios and the enormous data volumes.
In addition, the layout of the passive seismic monitoring network is key for the quality of data acquisition. We use the phase-1 of the dense array as an example to study the influence of the interstation spacing on the minimal magnitude of processed events. The  stations of our phase-1 dense array are roughly deployed within a rectangle of 30 km by 20 km (Figure 13). We select 25, 50 and 100 stations which are roughly in a uniform distribution out of the 187 stations to construct three sparser seismic networks covering the same area. We assume an earthquake can be reliably located and its event magnitude can be further analyzed if 8 P-wave and 8 Swave arrivals (P-and S-arrivals not necessarily from the same stations) are observed by a sparse network. With the increase in the station number, 21%, 60%, and 91% of total events (3961) monitored by the original phase-1 dense network of 187 stations are effectively monitored by the sparser networks with 25, 50 and 100 stations, respectively. Even with a relatively sparse network with 25 stations, the is 0.6, much lower than that of the CNNC. When the number of stations exceeds 50, the is reduced to , the same as the of the original dense array. However, it can be found that some seismic events greater than are omitted in the catalog of the network consisting of 50 stations, which means the estimated using the MAXC is biased (Figure 14). This is because the value using the MAXC method is determined by the magnitude bins with the highest values in the non-cumulative FDM (Figure 14b), which in the case of the sparser array with 50 stations are still around magnitude 0.0. The Goodness-of-Fit (GFT) technique (Wiemer and Wyss, 2000) is also used to reappraise the of the catalog recorded by the 50 stations, and the estimated  those by the original dense array, we can find the performance of these two networks is rather similar, with less than 10% of events smaller than omitted by this sparser network. Thus, the number of stations in the dense network can be halved in practice without seriously compromising the monitoring performance in terms of event magnitude.
As for the accuracy of magnitudes determined by the four different networks, we select the 812 common events recorded by the networks consisting of a varying number of stations, and calculate their magnitudes using Equation (16). Figure 15 shows the discrepancies in the event magnitudes between the sparser networks and the original dense network. In general, the accuracy of the event magnitudes determined by sparser networks is still quite satisfactory, with discrepancies no larger than 0.2. Also, the discrepancy decreases with increasing number of stations in the sparse network.

Conclusion
By analyzing the earthquake magnitudes monitored by a twophased dense seismic network covering the Changning−Zhaotong shale gas field, we show that the stations less than 10 km in epicentral distance tend to overestimate the local magnitude, and  more distant stations may underestimate the magnitude if the current national standard is used. This is because the amplitude calibration function does not properly reflect the rapid variation in amplitude for stations in close proximity to the events. Therefore, the current calibration function may not be suitable for dense seismic networks which have been or will be deployed to monitor the unconventional shale gas fields in the southern Sichuan Basin for the frequent seismic activities. To derive a proper and accurate formula for calculating the local magnitudes of earthquakes supposedly induced by hydraulic fracturing in the Changning− Zhaotong shale gas field, we analyze some 1,000 common events monitored by both CSN and our dense array, and determine the optimal fitting coefficients. The new local magnitude formula should be appropriate for the study region for stations with epicentral distance ranging from a few hundred meters to about 30 km, and is consistent with the current national standard. Using the calibrated formula, we determine the local magnitudes of 7,543 microearthquakes monitored by our dense network from February 28th to May 6th, 2019, processed manually. The minimal magnitude for the monitored events is −1.0, while the magnitude of completeness is further reduced to −0.1 in this area compared to 1.1 provided in the CNNC. In addition, we also use several machine-learning techniques to process the continuous waveform data, and further increase the number of detected and processed events from 7,543 to about 14,000. The in the region is even reduced to −0.3. The influence of station density on the minimal magnitude of the detected events is also investigated.
We believe many more studies using dense seismic networks will be conducted in the Changning shale gas field in the near future to better characterize the increasing seismicity and understand the correlation between hydraulic fracturing and fault reactivation, as well as the associated seismic risks and potentials. Therefore, our study will provide an important framework for calculating local magnitudes for dense seismic networks.

Data and Resources
The CNNC catalog is provided by China Earthquake Networks Center, National Earthquake Data Center (http://data.earthquake. cn). All the codes for magnitude calibration and calculation related to this study are accessible at the following website: https://zenodo.org/record/4660173#.YGc3yHUzakA.
To evaluate the individual contributions from the near, intermediate and far-field terms of a point dislocation source, we have performed a synthetic test in which stations are deployed around an earthquake of an arbitrary magnitude (we are only concerned with the relative amplitudes of the near, intermediate and farfields), with epicentral distances ranging from 1 km to 30 km at a 1 km interval at different azimuthal angles, as shown in the follow- are calculated by Equation (16), while are calculated by (Havskov and Ottemoller, 2010). ing Figure B1. We calculate the velocity amplitudes of the near, intermediate and far-fields observed by stations at different distances according to Eq. (4.29) in Aki and Richards (2002), averaged over four different azimuthal angles. The following figure

Earth and Planetary
shows that the near-field has negligible contribution compared to the intermediate and far-fields, and even the contribution from the intermediate field is much smaller than that from the far field, which decays as 1/r.  Figure B1. The variation of velocity amplitudes observed by stations at different epicentral distances. The inset shows the observation system, and the source is located at 1.5 km in depth. Stations are deployed at four different azimuthal angles, with a 1 km radial interval.