Bulk Crustal Properties and Layered Velocity Structure in Fujian, SE China: Constraints From P and S Receiver Functions

Multiple tectonic events since the Neoproterozoic era have framed the present-day lithosphere in the Fujian province affiliated with the eastern part of the South China Block. Comprehensive information of the crustal structure and bulk properties can aid to understand the geological features and tectonic processes of still much debate in this region. An attempt is made in this study to explore crustal thickness and internal velocities across Fujian using the teleseismic receiver functions (RFs). The H-V stacking of joint P and S RFs improves to simultaneously estimate crustal thickness, average Vp and Vs, and derived Vp/Vs ratio and bulk sound speed in three backazimuth sectors for each of 17 stations. Furthermore, a Neighborhood Algorithm nonlinear inversion of P RFs is employed to determine the layered structures of Vs and Vp/Vs beneath all the stations. Results indicate the crustal thickness varies from at most ∼35 km in northwest Fujian to 30–35 km in the inland mountains and 27–30 km in the southeastern coasts. The inferred Moho geometry is nonplanar or inclined across the Zhenghe-Dapu (ZD) and Changle-Zhaoan (CZ) fault zones, especially in the southern ZD fault area. The average Vp/Vs suggests that the crust is predominantly felsic in the Wuyi-Yunkai orogen and intermediate-to-mafic in the Cretaceous magmatic and metamorphic zones. A high-velocity upper crust along the coastline is revealed, which attributes to the Pingtan-Dongshan metamorphic belt. At the sites near the ZD fault zone, the intracrustal negative discontinuity occurs at a shallower depth of ∼15 km marking an abrupt Vs decrease into the low-velocity mid-to-lower crustal layer, probably linked to the closed paleo-rift basin remnants. The lower crust across the Fujian is generally characterized by relatively lower Vs and higher Vp/Vs (1.80–1.84) consistent with those of the mafic-ultramafic rocks, which do not support the proposed extensive magmatic underplating in the Late Mesozoic.


INTRODUCTION
The Fujian province on the southeastern coast of China resides on the Cathaysia Block, the southeast part of the South China Craton formed by collision with the Yangtze Block in the northwest during the Neoproterozoic era ( Figure 1A). Since the Paleozoic era, it had undergone several active tectonic events, such as the intracontinental orogenesis and collapse, westward subduction of the paleo-Pacific plate, lithospheric thinning, and extensive magmatism (Hsü et al., 1990;Li, 2000;Wang et al., 2013;Lin et al., 2021). The tectonic environment had thus changed from intraplate compressional to extensional stress regime in the Mesozoic (Li and Li, 2007) and then been under extension again in the Cenozoic (Qiu et al., 2018). These tectonic events had led to the formation of numerous uplifted and depressed blocks separated by the NE-SW and NW-SE oriented faults or sutures (Wang et al., 1993; Figure 1A), among which the major ones are the NNE-NE trending Shaowu-Heyuan (SH), Zhenghe-Dapu (ZD), Changle-Zhaoan (CZ) fault zones, and offshore Fujian littoral (FL) fault. The Wuyi-Yunkai orogenic belt located in the western Fujian is bounded by the SH and ZD suture zones while its eastern neighboring terrane is the Cretaceous magmatic downfaulted zone sided by the ZD and CZ fault zones (see Figure 1B). To the east of the CZ fault zone lies the Mesozoic Pingtan-Dongshan metamorphic belt (PDMB) along the coastline (Chen et al., 2002), which is delimited by the active FL fault in the western Taiwan Strait Zhang et al., 2020a). The FL fault marks the western boundary of the strongly rifted Cenozoic basins in the Eurasian continental margin that has been flexed down by the load of the Taiwan orogen to the east (Lin et al., 2003;Zhan et al., 2004).
The aforementioned geological blocks are characterized by the widespread Triassic to Cretaceous igneous rocks resulting from extensive magmatism (Zhou et al., 2006;Guo et al., 2012). The outcrop of the western Fujian mainly consists of Triassic and Jurassic granitic rocks while that of the eastern Fujian is comprised of Cretaceous granitic and volcanic rocks. The geological and geophysical studies have put forward several possible subsurface features that may prevail in the presentday lithosphere beneath Fujian, for example, replenished thin lithosphere (Ye et al., 2014;Zhang et al., 2020b), upper mantle lateral heterogeneity (Huang et al., 2010;Cai et al., 2015), deep crustal melting and mantle-derived magmatic underplating (Zhou and Li, 2000;Huang et al., 2011), intracrustal The upper right inset map displays the Yangtze and Cathaysia Blocks amalgamated to form the South China Craton. The 14 FEA and three BATS stations, located within our study region bounded by the cyan box, are denoted by blue and magenta triangles, respectively. Red lines delineate the faults mapped in the land (solid) and offshore areas (dashed) (Ding et al., 1999), while black lines draw the border of the Fujian province. A solid yellow dot near the Min River marks the location of Fuzhou basin. The NE-SW trending suture zones in the Fujian region are abbreviated: SH, Shaowu-Heyuan fault zone; ZD, Zhenghe-Dapu fault zone; CZ, Changle-Zhaoan fault zone; FL, the Fujian littoral fault. (B) The numbers of P and S RFs from individual backazimuth sectors (NW, NE, and SE) used for H-V stacking. Gray-colored sectors with smaller and larger radii indicate the respective numbers of P and S RFs in the three backazimuth ranges. The upper right inset map shows the epicentral distance and azimuthal distribution of teleseismic events providing usable P (solid black circle) and S (green open circle) RFs.
low-velocity zone (Cai et al., 2015;Zhang et al., 2018), lithosphere delamination with foundering of lower crust (Wang et al., 2013;Lin et al., 2021), crustal thinning and modification (Dong et al., 2020). However, more promising in-depth evidence for the nature and properties of these features is a prerequisite and worth being explored for better understanding of the geological evolution of the crust and upper mantle of Fujian and the earthquake-prone fault mechanisms. The receiver function based investigations for the Fujian region are mostly centered on the overall regional variation of crustal thickness and Vp/Vs ratios (Ai et al., 2007;Li et al., 2013;He et al., 2014), despite a few on the lithosphere thickness and upper mantle discontinuities (Ye et al., 2014;Huang et al., 2015), but they lack in seeking out the information of detailed crust and upper mantle structure.
In this study, we first employ P and S receiver functions (RFs) to explore both the azimuthal and lateral variation of the crustal thickness and bulk properties beneath available seismic broadband stations distributed across the Fujian province using the H-V stacking approach which was first introduced in detail by Rychert and Harmon, 2016 and later applied scrupulously to the Taiwan orogen in our previous study (Goyal and Hung, 2021). We next conduct a comprehensive nonlinear inversion of the stacked P RFs for each site to scrutinize the underlying shear velocity (Vs) structure of the crust. The results reveal the crust-mantle interface or Moho geometry derived from the RFs observed at three available backazimuth swaths, the pronounced intracrustal discontinuities and lowvelocity layers, and the depth-varying Vs and Vp/Vs ratio of the crust beneath each station. Lastly, we discuss their possible implications for the internal structures and tectonic processes of the southeastern China.

DATA AND RECEIVER FUNCTION PROCESSING
A total of 17 permanent broadband stations providing accessible waveform data are utilized to jointly evaluate the crustal properties and seismic velocity structure beneath the Fujian province (see Figure 1A). Most of these stations started operating in mid-2011, and since then, they have been maintained by the Fujian Earthquake Agency (FEA), while the remaining three located at nearshore islands belong to Broadband Array in Taiwan for Seismology (BATS) of the Institute of Earth Sciences (IES), Academia Sinica. We target seismic waveform records generated from the teleseismic earthquakes occurring between 2012 and 2021 with moment magnitude (Mw) greater than or equal to 5.7 and within the epicentral distance range of 30°-90°from the study region. The depths of events are limited to be less than 300 km to avoid the contamination of strong multiple surface-reflected P wave energy from deep-focus earthquakes (Wilson et al., 2006). To extract the P receiver function (RF) for each event-station pair, we first window the three-component seismograms starting from 30 s before and ending 120 s after the predicted P arrival and rotate them to the R(radial)-T(tangential)-Z(vertical) coordinate system. A zero-phase Butterworth bandpass filter of 0.05-1.00 Hz is applied to all the traces, following which the R-component trace is deconvolved from its Z-component counterpart using an iterative time-domain deconvolution scheme (Ligorría and Ammon, 1999). A Gaussian width (Gw) of 2.5 for the H-V stacking and 4.0 for the nonlinear RF inversion is chosen which yields the pulses with a dominant period of 1.6 and 1.0 s, respectively. The obtained RF is accepted for further analysis only if the variance reduction between the corresponding observed and predicted R-component waveforms is larger than 70%, thereby leading to nearly 330 events that have produced good-quality P RFs.
The S RF that comprises the S-to-p conversion from the Moho and its crustal multiples arriving before and after the direct S, respectively, is most prominent in optimal epicentral distances of 55°-85° (Yuan et al., 2006). Hence, we first window 200-s long waveforms centered at the predicted S arrival and rotate them from the Z-N-E to P-SV-SH coordinate system using the free surface transformation (Bostock, 1998). The direct S phases are identified and marked on the SV component using the AIMBAT software, which provides an interactive interface to measure phase arrival times from waveform cross-correlation (Lou et al., 2013). For the upkeep of quality, only those S waveforms that have signal-to-noise ratios (SNR) greater than three are retained. We utilize the extended-time multi-taper frequency-domain cross-correlation estimation (Helffrich, 2006) for deconvolution of the P from the SV component, and the deconvolved traces are then bandpassed between 0.03 and 0.50 Hz using a fourth-order Butterworth filter. We carefully inspected every calculated S RF to ensure that the direct converted and reverberated phase arrivals are coherent and have discernible amplitudes and correct moveouts varying with epicentral distance. Such selection procedure yields approximately a total of 100 S RFs per station. The azimuthal and distance distribution of the teleseismic events providing P and S RFs for our study are shown in Figure 1B.

Fundamentals of Methodology
The P RFs are fundamentally bedecked with direct P-to-s conversion (Ps) and associated reverberated phases (PpPs and PpSs+PsPs) due to the confrontation of an incoming P-wave with a lateral seismic discontinuity at depth. From a ray-theoretical point of view, the arrival times of these phases (t Ps , t PpPs , and t PpSs+PsPs ) relative to that of the direct P can be formulated as follows (Zhu and Kanamori, 2000): where H is the depth of the accountable seismic discontinuity, p is the horizontal slowness of an incoming teleseismic P wave, and V p and V s are the average compressional and shear wave velocity of the aggregate structure between the discontinuity and recording site, respectively. The above formulation presumes a high-frequency wave to encounter a flat discontinuity with sharp impedance or velocity contrast along with the laterally smoothly varying isotropic velocity structure. The crust-mantle boundary or the Moho is in general known as one of the most significant discontinuity featuring highamplitude converted and reverberated phases in P RFs. The above equations have been widely used to evaluate the crustal thickness (H) and V p /V s ratio or κ by means of the H-κ stacking technique of Zhu and Kanamori (2000), which performs a grid search for the optimal H and V p /V s ratio, given a priori assumed constant Vp of the crust, to reach the maximum value of a function defined as the weighted sum of amplitudes of the Moho related phases marked at the corresponding predicted arrival times relative to direct P. In reality, the Moho discontinuity can be inclined or nonplanar and behave as a gradual velocity transition (Endrun et al., 2005;Tang et al., 2011;Saikia et al., 2016). Moreover, the continental crust is usually multi-layered and anisotropic and may consist of a thick, low-velocity sedimentary cover on the top (Schulte-Pelkum and Mahan, 2014; Yu et al., 2015). Consequently, the P RFs may not exhibit coherent and prominent Moho converted and reverberated phases, which restrain the H-κ technique from determining a well-constrained global maximum solution of H and V p /V s .
On the other hand, the S RFs have been proved useful as the conversion phase (Sp) arrives earlier and the associated reverberations (SsPp, SsSp) later than the direct S, hence mitigating some of the stated difficulties while using the P RFs only. The equations for the arrival times of these phases (t Sp , t SsPp , and t SsSp ) are expressed as (Rychert and Harmon, 2016), It is important to note that the three arrival-time equations for P or S RFs could not uniquely determine all the three parameters (H, V p , and V s ) and, in addition, the secondary multiples, PpSs+PsPs and SsSp phases, are often less coherent and weakly ambiguous. Thus, only two crustal parameters (H and V p /V s ) are evaluated in the H-κ stacking of P or S RFs using the above equations. Because the other direct conversions and primary multiples (Ps, PpPs, Sp, and SpPp) in P and S RFs are usually more prominently identified, their arrival-time equations can be jointly incorporated into simultaneously robust estimation of all the three parameters in H-V p -V s space, which was first introduced as H-V stacking method (Rychert and Harmon, 2016) and later modified for mapping the complex Moho and crustal elastic properties beneath the Taiwan orogen (Goyal and Hung, 2021). In this H-V method, no a priori information for Vp is required, and the trade-off between H and V p /V s largely involved in the conventional H-κ analysis is much reduced (Goyal and Hung, 2021). The 3-D grid search for the first-order crustal parameters (H, V p , and V s ) can be carried out by maximizing the value of function (F) defined as where a is the amplitude at the respective phase arrival time, and w 1,.., 6 represents the weights of the corresponding phases.

Practical Approach
We conduct the joint H-V stacking of P and S RFs for each of the three main backazimuth sectors of teleseismic events separately at every station site. Since the obtained P and S RFs span over a wide range of backazimuths and epicentral distances ( Figure 1B), we categorize the usable RFs into three directional subsets: northwest (NW), northeast (NE), and southeast (SE) with their respective backazimuth intervals of 280°-320°, 20°-50°, and 110°-140°. Except for the southwest from which the usable RFs are considerably fewer and insufficient enough to carry out the H-V stacking protocol, we stack both the P and S RFs from each of the rest three backazimuth quadrants in 5°-wide bins of epicentral distance. For each station, a 3-D grid search is performed for individual backazimuth quadrants over a wide range of H, V p , and V s which are set to vary from 20 to 60 km, 5.5 to 7.0 km/s, and 3.0 to 4.0 km/s, respectively. In practice, we stack all the P RFs as one set, the S RFs as the other, normalize each set by the respective number of the stacked traces, and then seek the optimal combination of H, V p , and V s to maximum the function F defined in Equation 7. Because the secondary multiples, PpSs+PsPs and SsSp, are less coherent and difficult to discern for most of the stations, we choose the weighting scheme similar to that in Rychert and Harmon, 2016 as: w 1 0.25, w 2 0.20, w 3 0.00, w 4 0.30, w 5 0.25, and w 6 0.00 in the Equation 7. As the S RFs provide key constraints on Vp as demonstrated in Equation 5, we assign slightly higher weights for the S RF phases. In the grid search outcomes, we find the calculated F values in the 3-D H-V space are not sharply centralized around the global maximum value such that the volume for the 95th percentile isovalue contour is quite large for most of the stations, possibly due to the broader pulses in lower-frequency S RFs. This suggests that the solution corresponding to the maximum F value may not represent the most reliable and unbiased one and differs substantially that taken from the mean of statistically acceptable solutions around the maximum, as also reported in Rychert and Harmon, 2016. Therefore, instead of directly picking up the maximum value from the grid search, we determine the optimal solution of H, V p , and V s by taking the mean of the solution grids with F values over 95% confidence level and the associated uncertainties from their respective standard deviations. Besides, we also test other slightly different combinations of the weights for the four retained phases. They all yield similarly broad F-value distributions and comparable mean solutions within the 95% confidence regions, demonstrating that our H-V results are not affected substantially by the chosen weighting scheme.
To illustrate our H-V stack approach, we plot the contours for 90% and 95% confidence regions of the F values around the resulting global maximum projected onto H-V p , H-V s , and V p -V s planes using RFs from the SE quadrant for station DSXP as shown in Figures   More details of our H-V grid search procedure can be referred to Goyal and Hung (2021) in which we effectively applied it for determining average crustal Vp, Vs, and Moho depth beneath stations in the highly complex tectonic setting and orogenic structure of Taiwan. The results are remarkably consistent with those from the previous P RF analysis at several common stations and the depth-averaged Vp, Vs, and Vp/Vs over the vertical crustal column beneath each site calculated from established local tomography models. Moreover, a synthetic test was presented the supplementary figure in the published paper of Goyal and Hung (2021) to demonstrate that the H-V approach works capably for a site with a thick low-velocity sedimentary cover over a sharp intracrustal discontinuity that mimics the intra-continental crustal structure observed in Taiwan and Fujian. Here, we compare our findings of Vs, Vp, and H for a given example station DSXP by juxtaposing them with those from the two previous studies in South China based on the H-κ analysis of P RFs only  and the joint inversion of P RF and gravity data (Guo et al., 2019) in Figures 2C-E. Their results comply substantially with ours within the uncertainties except for the average Vp which was therein assumed to be fixed at a constant value of 6.3 km/s. All of these suggest that the employed H-V method is quite suitable for deciphering the first-order crustal features at sites beneath the Fujian province.

Inversion Modeling of P RF Stacks
The RFs are primarily sensitive to velocity contrasts across lateral discontinuities, although some investigations claim their potential to resolve the absolute velocities as well (Ammon, 1991;Svenningsen and Jacobsen, 2007). We adopt a nonlinear and derivative-free direct search method, named as the Neighborhood Algorithm (NA; Sambridge, 1999) to invert for the layered shear velocity (Vs) structure beneath each station from the P RF stacks. In this well-established and tested inversion scheme, the earth model is composed of six horizontal layers from the top, corresponding to sediment, basement, upper, middle, lower crust, and upper mantle, and each of them is characterized by four parameters: the layer thickness, Vs at the top and bottom depth of the layer, and constant Vp/Vs ratio in the layer, which totally makes up a 24-dimensional parameter space. The model is estimated to minimize the objective function defined by the chisquare misfit (χ 2 ) between observed and predicted data as follows: where N is the total number of data points in a stacked RF trace used for the inversion, O i and P i are the ith observed and predicted value of the RF, respectively, and σ i is the standard deviation calculated during the stack of individual RFs.
To resolve the layered structures beneath each station with higher resolution by the NA inversion, first of all, we select and stack the best quality, shorter-period P RFs which are calculated with Gw of 4.0 (equivalent to 1-s dominant period) and unambiguously consistent in a narrow range of epicentral distance and backazimuth. Thus, a total of five to nine stacked P RFs per station is obtained for a variety of distinct azimuth and distance bins, and each of them is then processed for the Vs and Vp/Vs models using the NA inversion approach. Each inversion run involves 800 iterations generating a total of 64,100 candidate velocity models and forward calculating the corresponding predicted synthetic RF and the objective function to be minimized, which measures L 2 -norm misfit between the observed and synthetic RF data. The NA inversion procedure is exemplified in Figures 3A,B for station DHTZ, showing the comparison between an observed stacked P RF and its synthetic counterpart predicted by the model at the minimum objective function ( Figure 3A), all the trial models (in grey) along with the 1000 best models (colored) and the best-fit model (thick red line) ( Figure 3B), and the monotonically decreasing χ 2 misfit with the increasing iterations that ultimately reaches a stable minimum ( Figure 3C). Finally, a laterally averaged Vs and Vp/Vs model is derived by averaging all the best-fit models resulting from the inversion of individual RF stacks ( Figures 3D,E). Since the NA approach implicitly assumes a laterally smooth layered structure beneath the targeted site, we exclude station YDFS from the current inversion procedure as its H-V results show considerably large variation of crustal thickness and average velocities among the three backazimuth sectors and will be discussed further in the next section.

Crustal Thickness and Average Seismic Velocities
As discussed in our earlier study (Goyal and Hung, 2021), the conversion points of Ps and Sp phases from useful epicentral distance ranges are located apart and off from the station by a horizontal distance of 6-11 and 21-30 km, respectively, given that the Moho lies at a depth of 35 km. Yet, despite that, their associated strong off-path sensitivity due to the inherently broad finite-frequency kernels (on the order of 15 and 35 km in crosspath width for Ps and Sp at a dominant period of 1.6 and 4 s, respectively) does partially overlap each other between the conversion points (Hung et al., 2000;Hung et al., 2001). Hence, the joint H-V stack of P and S RFs essentially yields the crustal thickness (H) and velocities (Vp and Vs) that are averaged laterally over the overlapping area about a dozen of kilometer wide toward the backazimuth quadrant. We thus employ this method to explore the variation of the bulk crustal properties over three different backazimuth sectors (280°-320°, 20°-50°, and 110°-140°) surrounding each station as shown in Figures 4, 5. Table 1 provides the location of each used station and obtained H-V results, including the determined thickness (H) and average Vp and Vs of the crust and the derived Vp/Vs ratio and Vb (bulk sound speed) with their associated uncertainties.
The H-V stack reveals that the crustal thickness varies gently from 27 km in the southeastern coast to 35 km in the northwestern end of the Fujian province. In western Fujian affiliated with the early Paleozoic Wuyi-Yunkai orogenic belt in South China, the crust has the thickness of 30-35 km and the low to normal Vp/Vs ratio of ∼1.66-1.77 typical for an intraplate continental crust with dominantly felsic composition. The thickest crust (34-35 km) is found in the northernmost Fujian with moderate Vp/Vs but higher Vp and Vs. The downfaulted magmatic belt in central Fujian shows the similarly thick crust (30-34 km), but having an increased Vp/Vs of 1.75-1.80 likely due to a decrement in Vs. Along the Fujian coastline, the crust exhibits the most undulating thicknesses varying from 27-30 km in the south to 29-33 km in the north and the overall slightly higher Vp/Vs than normal, except for the unusually high Vp/Vs (∼1.86) and much lower Vs (∼3.3 km/s) observed in few directional quadrants. Unlike greater variability in Vp values found in inland stations, those in the coastal stations are much close to 6.3 km/s. Some of the stations located near the Zhenghe-Dhapu (ZD) and Changle-Zhaoan (CZ) fault zones reveal quite different crustal thickness or Moho depths from the three backazimuth bins, implying a localized nonplanar Moho geometry. In general, the Moho is found deeper in the east of the ZD and the west of the northern CZ fault zone by about 2-4 km. The largest backazimuthal difference in the Moho depth over 4 km is obtained for station YDFS in the southern ZD. At this site, the crustal velocities estimated in the three different quadrants are contrasting as well, indicating the possibly distorted or disrupted Moho interface and substantial velocity contrast across the fault zone. The other noticeable difference in the Moho depth is observed at station MHZQ, where the crust in the Fuzhou Basin sampled by the RFs from the SE is shallower by 2.5-3 km than that from the NE and NW quadrants. Besides, for the two stations at the nearshore islands, VWUC in the eastern proximity of the littoral fault (FL) and KMNB just off the southeast coast of China, both exhibit a shallowing Moho toward the NE.

Thickness, Vs, and Vp/Vs of Layered Crustal Structure
The Fujian region has undergone several stages of tectonic activities that gave rise to the compositionally differentiated crustal strata cut by widespread faulted sutures. The layered intracrustal structures and properties can be hardly known from the average seismic velocities beneath the station sites determined by the H-V stacking. In order to gain a more comprehensive understanding of the crustal structure and its implications for regional geotectonics, the detailed velocity structure of the crust and uppermost mantle is required to be investigated. Here, we derive the Vs and Vp/Vs models beneath each station (except YDFS) from the stacked P RFs with Gw of 4.0 using the NA inversion approach. The inversion is performed individually for each RF trace calculated from the linear stack of carefully selected high-quality P RFs within a narrow range of backazimuth and epicentral distance. As the Ps piercing points at the intracrustal and Moho discontinuities are very close, less than 10 km to a station site and the sensitivity zone of high-frequency Ps phases dominated at 1 s is quite slim about 10 km wide, we integrate and average all the best-fit models inverted from each of the RF stacks over various backazimuth and distance bins to yield a robust estimate of the layered velocity structure beneath the station (e.g., Borah et al., 2016). Figure 6 and Table 2 shows all the obtained Vs and Vp/Vs models in various geological blocks along with the identified sediments colored in orange, the Moho, intracrustal positive (PD), and negative (ND) discontinuities marked by arrows. It may be noted that the most resolvable parameter in the RF inversion is the depth and sharpness (velocity contrast) of seismic discontinuities and then the absolute shear velocities of the crustal layers and their corresponding Vp/Vs variation.
The inverted shear velocity (Vs) models of the crust beneath the stations across Fujian, in general, resolve a PD followed by a ND as going deeper which grants us to identify four crustal layers: upper, middle, and lower crust separated well by the PD and ND, and the sediments, if any. The sedimentary layer is absent in most of the sites in the inland mountainous Fujian except at stations SXFK and MHZQ underlain by unusual, about 2-km thick sediments. In contrast, the crust beneath all the stations situated along the coast is overlain by 0.5-1-km thick sediments with Vs between 2.6 and 2.9 km/s ( Figure 6). The upper crust is generally 8-11 km thick with Vs ranging from 3.2 FIGURE 4 | Variation of crustal thickness between the stations estimated by the H-V stacking over the background Bouguer gravity anomaly map. The thickness resulting from the three backzimuth quadrants is shown by colored sectors with the radius also proportional to the corresponding thickness.
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 743163 to 3.4 km/s in the Wuyi-Yunkai orogen and Magmatic zone, and from 3.4 to 3.6 km/s in the coastal regions and all with a normal Vp/Vs ratio. We also notice an unusually thicker upper crust of ∼15 km beneath stations DSXP and XPSS in northernmost and southernmost coast of Fujian, respectively. Among all the models, those of the five stations (SXFK, YXBM, ZPXH, DHTZ, and MHZQ) situated in the inland central Fujian region around the ZD fault zone comprise of a uniquely shallow ND at the depth of 15-16 km ( Figure 6). The middle crust beneath these sites is only about ∼7 km thick, about one-half to two-thirds of the thickness (10-13 km) in the rest of sites. The seismic velocity property of the middle crust commonly pertains to a Vs of 3.6-3.7 km/s and a low-to-normal Vp/Vs. The lower crust throughout the Fujian province is found to be as slow Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 743163 9 TABLE 1 | Summary of crustal thickness (H) and average P-and S-wave velocity (Vp and Vs) of the crust including their estimated uncertainties resulting from the H-V stacking for all the three backazimuth sectors at 17 broadband stations.  Figure 7. The results obtained using the two independent methods are generally in good agreement with each other. Despite of some differences observed such as the average Vs at three stations, PCNP, SXFK, and YXBM, and Vp/ Vs at NPDK situated in western Fujian, they still match decently well within our estimated uncertainty range.

Variation of Crustal Thickness and Properties in the Fujian Province
The decrease in crustal thickness from 30-35 km in the inland Fujian to 27-30 km in the southern coast which lies the Pingtan-Dongshan metamorphic belt (PDMB) is readily observed from both the H-V stacking and NA inversion results. The findings are in good agreement with the previous RF studies mainly based on H-κ stacking of P RFs (Ai et al., 2007;He et al., 2014;Li et al., 2013), and largely consistent with the lateral variation of the Moho interface mapped from local earthquake travel-time tomography (Cai et al., 2015;Kuo et al., 2016) and joint inversion of P RF and gravity data (Guo et al., 2019). Moreover, the featured larger velocity jump across the Moho and higher uppermost mantle velocity in the south-central Fujian coast are in accordance with those seen in the tomographic results (Huang et al., 2010;Cai et al., 2015). Such crustal thinning from the inland to coastal region can be ascribed to the combined effect of the back-arc extension and lithospheric erosion due to the westward subduction and dehydration of the paleo-Pacific plate (Arcay et al., 2006;Niu et al., 2015). A consistently thinner crust (27-29 km) in the southern coastline than that in the northern (29-33 km) may demonstrate an important role of the NWstriking Min River fault in separation and interaction of the two blocks by regional stress adjustment as also proposed previously Ye et al., 2014). While considering the H-V stack results obtained in three different backazimuth sectors as a whole, a gently eastward and westward deepening Moho boundary exists across the Zhenghe-Dapu (ZD) and northern Changle-Zhaoan (CZ) suture zones, respectively, including a slight Moho uplift in the proximity of Fuzhou basin (Figure 4). These Moho undulations are plausibly a consequence of the Cretaceous magmatic belt bounded by the deep downfaulted ZD and CZ sutures. In addition, the Moho is found to be substantially nonplanar beneath the southern ZD fault (station YDFS) where the seismicity is most abundant and reaches deeper into the mid-lower crust, suggesting that the ZD fault is deeply rooted and highly active in this region (Lin et al., 2021). The significant variation of crustal properties and large historical earthquakes recorded in this region compels for a further in-depth investigation of fault mechanism and stress pattern that prevail in the crust.
The Vp/Vs is widely accepted to be more effective than Vp or Vs alone for inferring the changes in chemical composition, temperature, and pore fluid of rocks (e.g., Christensen, 1996). The resulting average crustal Vp/Vs ratio is lower in the Wuyi-Yunkai fold belt (1.68-1.73) and comparatively higher in the Cretaceous magmatic zone (1.74-1.80). Such discernible Vp/Vs variation may mostly reflect the different/distinct compositional characteristics of crustal rocks beneath these regions, suggesting that the crust in the Wuyi-Yunkai orogen is more felsic and that in the magmatic downfaulted belt is intermediate-to-mafic in nature. As the peraluminous granitoids and Mesozoic volcanic rocks are widespread in the Fujian province, they are chemically characterized by high silica content and calc-alkaline rich plagioclase feldspar, respectively (Zhou and Li, 2000). From the compilation of laboratory-based seismic velocity measurements of igneous rocks at the P-T conditions at crustal depths (Christensen, 1996), either a decrease of silica or an increase of calc-alkaline content could lower the Vp/Vs ratio significantly, which offers the most probable explanation for the increased Vp/Vs in the eastern Fujian.
In the coastal and nearshore region, though the average Vp/Vs ratios obtained from different backazimuth sectors mostly fall within the normal-to-high range of 1.73-1.79, a few stations like LYJJ and XPSS yield considerably higher values (∼1.83-1.87) which seem to appear intermittently and irregularly with no specific locations and backazimuths ( Figure 5A). This indicates the crust at the shore area is laterally inhomogeneous with the change in rock composition sporadically from being intermediate to highly mafic. Such diverse characteristics can be also noticed in the velocity inversion model, discussed later in the next section. The unusually increased Vp/Vs along with the decreased Vs at  few locations may also be related to a high geothermal gradient as there are many hot springs along the Fujian coastline. Not only is the Vp/Vs derived from the H-V stacking but also the bulk sound speed (Vb) being directly linked to the bulk modulus or incompressibility of rocks can also provide a measure of the crustal resistance to compressional or extensional tectonic stress. The average crustal Vb is found to be remarkably low in the Wuyi-Yunkai orogenic belt near the ZD fault zone and consistently higher in the magmatic and coastal zones. We thus speculate that the crust to the west of the ZD suture has become mechanically weaker in response to the subsequent collisional and extensional tectonics since Paleozoic; while it remains relatively strong, less deformed in the eastern Fujian area. Interestingly, the characteristics of the lowest crustal Vb and Vp/Vs in the ancient Wuyi-Yunkai orogen are similar to that observed in the young central Taiwan orogen with the most thickened, compressible crust (Goyal and Hung, 2021).

Layered Crustal Structure and Its Implications
A clear positive discontinuity (PD) corresponding to a sudden increase in shear velocity with depth pertains throughout the Fujian province mostly occurring at depths of 8-11 km, and we identify it as the interface between the upper and middle layer of the crust. The slower upper crustal layer atop the intracrustal PD in inland Fujian is similarly found in the joint inversion results of P RFs and surface-wave dispersion with Vp constraints (Deng et al., 2019) where it was explained by the formation of the granitoids, which would eventually be exposed on the surface after long-time erosion. On the other hand, the resulting negative velocity contrasts (ND) show different degrees of Vs reduction and variable depth range (15-24 km). In our view, the explicitly pronounced positive interface (PD) at 8-11 km depth range can be possibly related to the Conrad discontinuity commonly found in the continental regions and correspond to the transition from the upper crust dominated by more felsic granitoids to the midto-lower crust with more mafic rock composition. Along the coastline, between the CZ and FL fault zones, the shear velocities of the upper crust are consistently faster than those in the inland Fujian by ∼0.2 km/s. Such anomalously high-velocity feature in the upper crust has also been found in the tomography results (Cai et al., 2015;Kuo et al., 2016;Zhang et al., 2018), and it resides within the proximity of PDMB where the crustal rocks are mainly composed of granites and high-pressure amphibolite-facies metamorphic rocks (Chen et al., 2002;Liu et al., 2012). The relatively uniform, flat-layered sediments of ∼0.5-1 km thick are deposited along the coastline while they are seemingly absent in the inland Fujian, except for the 2-km thick sedimentary layer found at two sites (MHZQ and SXFK) which might be related to the adjacent Fuzhou basin and river valleys of the Wuyi-Yunkai mountains ( Figure 6). On the other hand, the shear velocities are nearly uniform in the middle crustal layer of 3.6-3.7 km/s along with the low Vp/Vs (1.68-1.74), possibly attributed to the less adulterated felsic rocks in the ancient orogenic belt. At the base of the middle crust, there exists a negative velocity gradient (ND) at the depth range of 20-24 km except for that beneath the five sites in the southcentral Fujian around the ZD fault zone is located at the shallower depths of 15-16 km. The low-velocity anomalies in the depths of the mid-to-lower crust have been reported in this region from the tomography studies as well (Cai et al., 2015;Zhang et al., 2018), though the 3-D tomographic images show they are distributed more irregularly at varying depths. Since RFs are mainly sensitive to the velocity contrast between layered structures, we detect that this slow zone appears from ∼15 km depth and extends down to the base of the crust, and has a normal to slightly higher Vp/Vs between 1.77 and 1.83. The closure of a paleo-oceanic rifted basin as a result of the Indosinian collision has been proposed from the observations of the symmetrical step-like velocity structure (Kuo et al., 2016) and the exposed surface rocks near the ZD fault zone such as ophiolite mélange, pelagic sediments, and basaltic debris (Li, 2013). Our low-velocity layer present in the mid-to-lower crust with high Vp/Vs complies well with that of the oceanic basin remnants (Hyndman, 1979). Intriguingly, a recent study based on the pre-stack migration of P RFs found a similar low-velocity FIGURE 7 | Comparison of our results of the Moho depth, average crustal Vs, and Vp/Vs derived from the inversion models and H-V stacking for each station, arranged according to its geological setting in the same order as that in Figure 6 and Table 2. The open black squares represent the results from the NA inversion while the red-, green-, and blue-filled circles show those from the H-V stacking method for the quadrants of NW, NE, and SE backazimuths, respectively. mid-crustal layer which is laterally broadly distributed over the southeastern China block and vertically varies from 10 to 16 km depths (Zhou et al., 2020). However, they deduced a lower Vp/Vs ratio (∼1.66) for this layer and thus related it to be trapped fluids released from dehydration of the subducted paleo-Pacific plate. Nevertheless, a conspicuous low-velocity layer in the mid-tolower crust prevails in the southeastern China margin near to the ZD suture zone and a more comprehensive seismic, and other geophysical evidence for this feature is yet to be explored in order to explain its origin.
Though the interpretation of the earth's gravitational field for the depth, shape, and density of causative bodies is theoretically nonunique, the Bouguer gravity anomaly corrected for the known terrain attraction has been commonly used to deduce the variation of crustal thickness and density distribution. As extracted from the global gravity model of WGM 2012 (Bonvalot et al., 2012), the Bouguer gravity map in the Fujian region reveals the lateral variation trending NE-SW roughly parallel to the coastline (Figure 4), showing regional lows in the inland Fujian and highs in the coastal and eastern offshore areas. Such variation most probably implies either the thinning of the crust or the increase of crustal density or both toward the coast which overall match the crustal features resulting from our RF study. Namely, compared to the crustal structure found in the coastal metamorphic belt, a thicker crust along with a lower Vs (usually lower density) upper crust is revealed in the inland Fujian. In addition, a thicker, low-Vs layer in the mid-to-lower crust around the ZD fault zone may also contribute to the lower Bouguer anomaly there.
There are a few sites along the coastline that exhibit a distinct Vs structure. For instance, stations XPSS and DSXP located far apart in the northernmost and southernmost Fujian coast, respectively, both similarly consist of a sharp velocity increase across the sediment-basement interface and thicker upper crust (∼15 km thick) with nearly constant, relatively higher Vs (∼3.5 km/s), while the others (except for PTMZ) yield a thinner upper crust (≤10 km thick) usually with a gradational increase of Vs with depth. Unlike the velocity models for the rest of stations comprising a low-velocity sediment cover and a distinct positive discontinuity between the upper and middle crust, that for station PTMZ in the central coast situated at the exposed surface basement rock reveals a transitional decrease in Vs from the lower upper crust to the middle crust and a broad low-velocity zone in the middle crust. These seismic signatures are likely associated with the very locally distributed volcanic and intrusive rocks formed from the frequent magmatic and extensional events in the Mesozoic (Liu et al., 2012;Wang et al., 2017).
Based on geochemical and petrological studies (Zhou and Li, 2000;Zhou et al., 2006;Jiang et al., 2009;Huang et al., 2011), the widespread granitoid rocks are largely believed to be the product of crustal material melting associated with the asthenosphere upwelling and thinning of lithosphere and/or magmatic underplating at the base of crust. The lower crustal layer is found to have low shear velocities of 3.5-3.7 km/s and higher Vp/Vs (1.80-1.84) throughout the Fujian region that deduces compressional velocities of 6.37-6.80 km/s. These velocities are not fast enough to represent those of the underplated rocks which typically exceed 4.0 km/s for Vs (Saikia et al., 2017) and 7.0 km/s for Vp (Thybo and Artemieva, 2013). Rather, they are more consistent with the lower crustal rocks containing a mafic-ultramafic composition. The magmatic underplating frequently induces remelting of the lower crust (Fountain, 1989;Ma et al., 2012); therefore, it could also potentially generate the extensive granitoid magmatism (Wei et al., 2020). Such remelting could have reduced the lower crustal velocities despite remaining the Vp/Vs approximately unchanged. Alternatively, the thinning of lithosphere accompanied by asthenosphere upwelling may also explain the late Mesozoic magmatism (Deng et al., 2019;Lin et al., 2021).

CONCLUSION
We perform a thorough analysis of P and S RFs by employing the H-V stacking and nonlinear inversion methods to investigate the thickness, average seismic properties, and layered velocity structure of the crust beneath the Fujian province in southeastern China. The resulting crustal thickness varies from 30 to 34 km in the inland Fujian, being thickest about 34-35 km in the northwestern end and Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 743163 gradually thinned to 27-30 km in southern coastal region. The Moho boundary is found to be nonplanar or inclined across the active Zhenghe-Dapu (ZD) and Changle-Zhaoan (CZ) fault zones, with the most distorted geometry across the southern ZD fault, implying the fault roots deeply into the crust. The average crustal Vp/Vs reflects a more felsic composition in the Wuyi-Yunkai orogen (1.68-1.74), while more close to the intermediate-to-mafic characters in the Mesozoic magmatic and metamorphic coastal zones (1.73-1.81). The inverted Vs and Vp/Vs models reveal two distinct intracrustal discontinuities with an abrupt positive and negative velocity change with depth throughout the Fujian region, which mostly occurs at approximate depths of 8-11 and 20-24 km, respectively. The models manifest a 0.5-1 km thick sediment layer underlain by a relatively fast upper crust in the metamorphosed belt along the coastline. Beneath the sites in close proximity to the ZD fault zone, the negative discontinuity rather appears shallower at ∼15 km depth as the top boundary of a thick mid-to-lower crustal low-velocity zone. The lower crustal layer across Fujian is observed with lower Vs (3.5-3.7 km/s) and higher Vp/Vs (1.80-1.84) in accordance with the seismic properties of mafic-ultramafic rocks. A comprehensive interpretation of our main findings stated above is summarized and illustrated through a 3-D schematic diagram presented in Figure 8.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/ restrictions: the waveform data used in this study are made available through the Fujian Earthquake Agency (FEA) in