CCREM: New Reference Earth Model From the Global Coda‐Correlation Wavefield

The existing reference Earth models have provided an excellent one‐dimensional representation of Earth's properties as a function of its radius and explained many seismic observations in a broad frequency band. However, some discrepancies still exist among these models near the first‐order discontinuities (e.g., the core‐mantle and the inner‐core boundaries) due to different data sets and approaches. As a new paradigm in global seismology, the analysis of coda‐correlation wavefield is fundamentally different from interpreting direct observations of seismic phases or free oscillations of the Earth. The correlation features exist in global correlograms due to the similarity of body waves reverberating through the Earth's interior. As such, there is a great potential to utilize the information stored in the coda‐correlation wavefield in constraining the Earth's internal structure. Here, we deploy the global earthquake‐coda correlation wavefield as an independent data source in the 15–50 s period interval to increase the Earth's radial structure constraints. We assemble a data set of multiple pronounced correlation features and fit both their travel times and waveforms by computing synthetic correlograms through a series of candidate models. Misfit measurements for correlation features are then computed to search for the best‐fitting model. The model that provides an optimal representation of the correlation features in the coda‐correlation wavefield is CCREM. It displays differences in radial seismic velocities, especially near the first‐order discontinuities, relative to previously proposed Earth‐reference models. This is the first application of the earthquake‐coda correlation wavefield in constraining the whole Earth's radial velocity structure.

inner-core boundary (ICB). The structures near these boundary layers can provide critical constraints on understanding the Earth's interior dynamics as well as heat and material transport throughout geological time.
Traditionally, reference Earth models were constructed based on either short-period body-wave observations or long-period normal-mode data. These two types of data could provide different constraints on Earth's physical properties. Therefore, reference models are usually constructed for specific practical use but show limitations for other purposes. By analyzing travel times of major body-wave phases, several travel-time reference models have been established, such as iasp91 (Kennett & Engdahl, 1991), sp6 (Morelli & Dziewoński, 1993), ak135 (Kennett et al., 1995), and ek137 (Kennett, 2020). The iasp91 model is proposed to yield more accurate earthquake locations, whereas sp6 is designed to be a closer representation of the globally averaged structure. The significant motivation behind developing both ak135 and ek137 models was to represent the core phases' travel times better. The models derived from body-wave travel times can effectively locate global earthquakes and are appropriate references for body-wave tomographic studies. However, apart from seismic velocity information, the density distribution in these models cannot be tightly constrained (Kennett, 2020).
In contrast, to resolve the density distribution for the whole Earth, an alternative approach is to use eigenfrequencies of normal modes with/without body-wave data (e.g., Dziewoński & Anderson, 1981;Gilbert & Dziewoński, 1975;Jordan & Anderson, 1974). Among these models, PREM (Dziewoński & Anderson, 1981) has been an efficient reference model for the last 40 years. Nevertheless, because of the low frequencies of the available modes (<0.0125 Hz), the resolution of structures in depth is limited (Kennett et al., 1995). Therefore, models derived from normal-mode data are commonly used for long-period studies.
All these proposed models are similar and represent many aspects of the seismic wavefield (Kennett, 2020). However, some discrepancies in P-wave velocity structures still exist among these models, especially near the internal boundaries ( Figure 1). A plausible reason is that these reference models are constructed based on different categories of data sets and methods, which have different sampling sensitivities and resolutions to the structures. The inconsistencies in velocity profiles have thus motivated many studies on refining the Earth's velocity profiles, especially in the lower mantle and core (e.g., Kaneshima, 2018;Ohtaki et al., 2012;Robson & Romanowicz, 2019;Ruff & Helmberger, 1982;Song & Helmberger, 1995;Souriau & Poupinet, 1991;Tanaka, 2007). However, no consensus on the velocity profile in these regions has been satisfactorily reached. Here, we analyze a data source from earthquake-coda correlation wavefield instead of using regular body-wave travel time or normal-mode data to resolve Earth's radial velocity structure discrepancies.
The earthquake-coda correlation wavefield is a mathematical expression of the seismic wavefield presented as a 2-D global coda cross-correlation stacks known as a global correlogram (for a review, see . The global correlogram contains a wealth of prominent and stable correlation features, which exhibit noticeable similarities with regular seismic wavefield in time-distance stacks (e.g., Boué et al., 2014;Phạm et al., 2018;Poli et al., 2017;Ruigrok et al., 2008). However, some puzzling features not present in regular seismic wavefield were also observed in the correlograms (e.g., Boué et al., 2013;Pham et al., 2018) and were dubbed "spurious." Poli et al. (2017) attributed the seismic-phase-like features to the interference of high-order normal modes. Pham et al. (2018) used the ray theory to explain both the features that resembled seismic phases and "spurious" features that had no adequate explanation, as the interaction of many pairs of phases with the same slowness at a receiver pair. Kennett and Phạm (2018) extended the formalism in the context of the generalized ray theory. The coda wavefield has been shown to be made of energy reverberating in the great-circle plane (Poli et al., 2017;Sens-Schönfelder et al., 2015). These features can be further dissected into different constituents attributed to variable cross-terms of reverberating body waves through the whole Earth (Poli et al., 2017;Wang & Tkalčić, 2020a). Therefore, the correlation features in the global correlogram can be regarded as fingerprints of Earth's interior. Analysis of these features could place tight constraints on the whole Earth's radial structure.
The correlation wavefield has emerged as a powerful technique in global observational seismology in the last several years. Based on early understanding, a series of deep body-wave phases were identified and used to interpret Earth structures (e.g., Huang et al., 2015;Poli et al., 2017;Wang et al., 2015;Wang & Song, 2018;Xia et al., 2016). Resting on a new understanding of how correlation features form, Tkalčić and Phạm (2018) detected shear waves in the Earth's inner core. In parallel, deeply sampling body waves have also been identified from microseisms (e.g., Li et al., 2020;Poli et al., 2012) and used to image deep Earth (Poli et al., 2015;Retailleau et al., 2020). Furthermore, a formation theory for earthquake coda-correlation has been confirmed by observing correlation-feature constituents, and feasibility for correlation tomography has been demonstrated (Wang & Tkalčić, 2020a, 2020b. Tkalčić and Phạm (2020) have recently shown that individual large earthquakes are sufficient in creating a high-quality global correlogram. Thus, numerical simulations of the correlation wavefield based on several high-quality individual events are computationally affordable. These advancements have shown promising potential and laid the foundations for using correlated body-wave signals to study the Earth's structures through the correlation wavefield .    Motivated by recent developments described above, we take a step forward in utilizing the coda-correlation wavefield in this study. As an independent data set, we intend to place constraints on the Earth's radial velocity structures based on coda-correlation wavefield observations in the 15-50 s period range, which is between the periods used in constructing PREM and ak135/ek137 models. First, the global correlogram is built from stacking cross-correlations of late-coda recordings from ten selected large earthquakes. Using correlation features in the correlogram as our observations, we search for the best-fitting models by comparing the synthetic and observed data. Then we compare our optimal model with a set of well-known models and demonstrate the implications of the optimal model for Earth's dynamics. This is the first attempt to constrain spherically symmetric Earth's velocity structure from an approach not based on direct observations of seismic phases or Earth's free oscillations.

Observations
Unlike the traveltime data used in previous studies, we choose a set of prominent correlation features in the global correlogram as our observations. Not only the travel time but also waveform information of these correlation features are utilized to constrain the Earth's radial structure.

Construction of Global Correlogram
First, the global correlogram is constructed using only ten selective large earthquakes (Mw ≥ 7.0) (Table S1) instead of a large number of earthquakes, as was the case in early studies. The ten earthquakes, showing either normal or reverse focal mechanisms with short-duration source time functions ( Figure S1), are chosen from the catalog of high-quality events presented in . With the efficient release of body waves along the Earth's radius, these individual events have been demonstrated to be sufficient in creating a high-quality global correlogram. Their summation results in a global correlogram with equal or better quality compared to the correlogram stacked over many events (e.g., as in Pham et al., 2018). Besides, with these ten carefully selected events, it becomes computationally feasible to generate a synthetic correlogram.
We select the vertical component recordings because (a) steeply reverberating waves with near-zero slowness are dominant in the late coda recordings while regular seismic phases with larger slowness quickly fade away with time after the origin time; (b) prominent features in global vertical-component correlograms are formed due to many cross-terms of near-vertically reverberating body waves. In particular, large events with reverse or normal focal mechanisms and simple source-time functions radiate energy steeply downwards, contributing to the most prominent cross-correlations of reverberating waves ; (c) many correlation features derived from the cross-correlation of vertical component seismograms contain both P-and S-wave propagation legs due to the energy partitioning between P and S waves at the internal boundaries.
After removing instrument responses, the seismograms are decimated at a 1 Hz sampling rate. We then cut out the continuous late-coda recordings (3-9 h after the event origin time) and process the data following the procedures described in Bensen et al. (2007), Phạm et al. (2018, which include temporal normalization and spectral whitening methods. Subsequently, cross-correlations between receiver pairs are computed for all globally available stations ( Figure 2a). For some events, the number of stations is over 2,000. To save computational times in the later simulation, we reduce the number of stations by choosing a single station in a 0.5*0.5° meshing element on the Earth's surface. In this way, the number of recordings for each event can be reduced by about 20%. The total number of receiver pairs relative to inter-receiver distance on a global scale is shown in Figure 2b. After applying a band-pass filter of 15-50 s, we then construct the global correlogram by linearly stacking all the filtered cross-correlations with a bin size of 1°. In Figure 2c, we present the global correlogram as a function of angular inter-receiver distances for the first two hours after the correlation-origin time. A wealth of prominent correlation features (e.g., Boué et al., 2013;Phạm et al., 2018;Poli et al., 2017) can be visually observed. The naming convention and abbreviations are detailed in Tkalčić and Phạm (2020). It is worth noting again that these features are not "reconstructed" body-wave Green's functions but have been demonstrated to correspond to many cross-terms of multiple reverberating body waves through the whole Earth with a common ray parameter (Kennett & Phạm, 2018;Phạm et al., 2018;Poli et al., 2017). These body-wave cross-terms can enhance illumination of the Earth's interior that is not well sampled via direct body-wave observations .
For each single-event correlogram of the ten selected events, we observe slight variations in the travel times and amplitudes of the correlation features due to 3D heterogeneities of the Earth. However, the process of summing correlograms effectively smooth out 3D heterogeneity effects within the Earth. On a more fundamental level, any given coda-correlation feature is generated by multiple body-wave cross-terms (constituents), which sample the Earth along fundamentally different paths (Wang & Tkalčić, 2020a). 3D heterogeneity effects are thus smoothed out due to stacking of all constituents sampling different Earth's volumes. Moreover, such effects are minimized by stacking many thousands of receiver pairs in different locations and binning them in inter-receiver distance bins, reducing Earth's ellipticity effects on the correlation features.

Correlogram Features Selection
We select most of the labeled features in the catalog assembled by Tkalčić and Phạm (2020). We complement them with some unlabeled, late-emerging features in the correlogram. We develop an approach to determine the time windows and distance ranges for each correlation feature ( Figure 3). First, we visually choose relatively broad time windows and distance ranges for individual correlation features ( Figure 3a). For each feature, the middle trace is selected as our initial reference trace. We then calculate the correlation coefficient (CC) between the reference trace and the neighboring traces for a 100-s window, including the feature signal. The time window length is set as twice the longest period (50 s   is larger than 0.8, we stack the two compared traces as our new reference trace and repeat the above process until the CC value does not meet the criterion. In this way, the initial distance range could be narrowed down to the range where all the feature signals are generally coherent and clearly expressed ( Figure 3b). The selected prominent features can be confirmed by the prominent energy displayed in the slowness-time domain using the phase-weighted stacking method (Schimmel & Paulssen, 1997) ( Figure S2). Additionally, we only keep the features that emerge in more than five traces after the above selection process. We then bound the prominent feature on each trace in the 100-s time window based on the slant stacking method (Davies et al., 1971;Rost & Thomas, 2002). In total, we select 71 prominent correlation features (blue lines in Figure 4) in the correlogram as the observed data.

Measurements of Waveform Fit
The most direct way to derive a spherically symmetric Earth model is to linearize the problem and invert all available data in a least-squares sense. However, we cannot carry out an inversion because the exact derivation of sensitivity kernels of correlation features for Earth's physical parameters is complex (e.g., Sager et al., 2018). Moreover, it is computationally expensive to simulate and post-process so long-lasting coda waveforms in the inversion. Therefore, we use a grid-search method to find the best-fitting model by comparing the synthetics with observed features through a series of candidate models.
To quantitatively express the fitness of each candidate model, we construct three measurements of fit. The CC, phase correlation coefficient (PCC), and L2-norm misfit values are computed between the observed and synthetic correlation feature signal for a particular model for each trace. PCC is used here as a complementary criterion for the goodness of fit because it is not an amplitude-biased measurement, and it keeps waveform coherence (Schimmel, 1999). Besides, these measurements can inherently account for the measures of time variations between observations and predictions. The averaged CC, PCC, and L2-norm misfit values are computed for each feature, and a measure of overall performance is provided by summing up these three values for all selected features, respectively. Large CC and PCC values and low L2-norm misfits indicate small time-variations and high waveform similarity between observed and synthetic correlation features. Then we use the total summed-up values as a fit criterion to search for the best-fitting models.

Model Construction
We construct the candidate models in two steps. The first step is building up all the candidate models using a weighted combination of four base models. These base models are PREM (Dziewoński & Anderson, 1981), ak135 (Kennett et al., 1995), PREM with wave speeds reduced by 5%, and ak135 with wave speeds increased by 5% (Figure 5a). The ak135 model and PREM are chosen as base models due to their wide use as reference models in the seismological community in the past few decades. Here we only perturb the P-, S-wave velocities and fix the density perturbations as zero since, empirically, the seismic velocity plays a dominant role over density in affecting the correlation features. For simplicity, we initially fix the CMB and ICB to PREM values in candidate models and rule out the discontinuities in the upper mantle considering the discontinuities in the upper mantle for PREM and ak135 vary by ∼10 km. We also fix the attenuation to PREM values.
For mathematical simplicity, we parameterize the velocity and density structure in the base models as piecewise cubic polynomials in radius. Then velocities or densities at given depths in each model are calculated as the weighted sum of all corresponding values in the base models. The weight for each base model ranges from 0.0 to 1.0 with a step interval of 0.2. One inherent condition is that the sum of the four weighting factors should be fixed as 1.0. Based on this grid-search approach, we can cover a wide range of different Earth's radial structures as our candidate models (gray lines in Figures 5b and 5c). In addition to these models, we also take into account some previous models in the simulation, such as PREM2 (Song & Helmberger, 1995), EPOC (Irving et al., 2018), and ek137 (Kennett, 2020). A complete calculation of summed-up CC, PCC, and L2-norm misfit values for all correlation features is performed for each model. Models with the largest CC, PCC, and smallest L2-norm values are chosen as the optimal models that best fit the coda correlation data. Note that all the unknown features are denoted by a character "xx" plus a number. Moreover, if the feature expands over a long-distance range (>50°) or has a cusp, we split the distance range into several parts. The feature is then named xx-1, xx-2, or xx-3 (xx represents the feature name).
In the second step, we further refine and modify the obtained optimal models to better fit the correlation features. These modifications include adding other discontinuities (i.e., the Moho and the upper mantle discontinuities), varying the ICB depths, and using a smaller interval step in the grid-search approach. Because we derive a reference model from the coda-correlation wavefield, we name the final optimal model CCREM (Coda Correlation Earth Reference Model).

The Optimal Model (CCREM)
To generate the synthetic correlogram, we first compute 9-h synthetic seismograms for the 10 selected events with the source-receiver configurations corresponding to the recorded data using Yspec (Al-Attar & Woodhouse, 2008). The moment-tensor solutions and source-time functions are obtained from the Global Centroid Moment Tensor catalog (Ekström et al., 2012) and SCARDEC catalog (Vallée & Douet, 2016). The synthetic late-coda waveforms are processed in the same way as the observations. We further normalize the time-windowed waveforms to exclude the attenuation effects, considering the significant uncertainties in Earth's attenuation structure.
Among all considered previous models (e.g., PREM, PREM2, ak135, EPOC, and ek137), our synthetic results show that the ek137 model provides the best fit for the selected correlation features. This is because ek137 is developed to improve fits for the travel times of outer-core sensitive phases (Kennett, 2020). Therefore, we choose the ek137 model as our reference model in the following sections to demonstrate that the CCREM provides an overall best fit for selected correlation features. Besides, our simulation results show that the PREM2 model outperforms PREM in fitting these correlation features since PREM2 offers a better match to different types of core phases in terms of PKP differential travel times, amplitude ratios, and waveforms (Song & Helmberger, 1995). Considering the latter improvement, some models constructed in Section 3.2 are further modified and tested by replacing the P-wave profile of PREM with that in PREM2.
MA AND TKALČIĆ 10.1029/2021JB022515 9 of 20 Based on the fit quantification described in Section 3.1, we display the best-fitting model in Figures 5b and 5c. Our optimal model, CCREM, has the largest summed-up CC, PCC values, and the smallest L2-norm misfits. In addition, to statistically estimate the significance of the CCREM, we evaluate the fit improvement to the data using the Wilcoxon signedrank test (Wilcoxon, 1945). Table 1 compares CCREM with three previously proposed reference models (PREM, ak135, ek137) via CC and PCC values. The test using CCREM as a reference model indicates that at a significance level of 5%, all the mean values of CC and PCC in PREM, ak135, and ek137 are smaller than those of CCREM. This implies that CCREM provides a significantly better fit to the correlation features than the other three models. The details on the test used to determine whether a model is significantly better than another model are in Text S1. Similar results can also be obtained by using the paired student's t-test (Table S2). Furthermore, we present frequency histograms of fits for the CCREM and ek137 models in Figure 6. Note. If the p-value is smaller than 0.05, we reject the hypothesis at a 5% significance level in support of the alternative hypothesis. Here, the parameter u represents the difference in CC/PCC values for all correlation features between the two models. The two hypotheses are: (a) u equals zero; (b) u is larger than zero, respectively. More details are described in Text S1. CC, Correlation coefficient; PCC, Phase correlation coefficient. selections on determining the best-fitting model. Each time, we calculate the CC, PCC, and L2-norm misfit values for the resampled features and obtain the optimal model with the largest CC, PCC, and smallest L2norm misfit values. Our bootstrap results show that out of 200 times, the CCREM stands as the best-fitting model in terms of CC and L2-norm values for more than 150 times. This means that the obtained best-fitting model (CCREM) is almost independent of the correlation feature selections. However, we note that another candidate model shows a slightly better performance than CCREM in PCC values. We do not choose this model as the best-fitting one because PCC is only used as a complementary criterion to CC and L2-norm misfit measurements.
On average, the CCREM CC value for each correlation feature is 0.711, and the PCC value is 0.532. A comparison between predictions from the CCREM and observations for two selected correlation features (PcP* and K4*) is shown in Figure 7. These features can be generally matched quite well in terms of both travel times and waveforms. To facilitate visualizing the waveform fit improvement of the correlation features, we further compare waveforms of these two features for CCREM and three reference models (PREM, ak135, and ek137) in Figure S3. We present the enlarged sections of the waveform comparisons for the four models in Figure 8. In addition, more waveform comparisons for other correlation features are displayed in Figure S4. We should keep in mind that the waveform fit is not improved for all the correlation features, and about 20% of the selected features are fit slightly worse for the CCREM (Figure 9). This is possibly due to different sampling sensitivities of various features to the radial structures. Despite reduced fitness for some MA AND TKALČIĆ 10.1029/2021JB022515 11 of 20 features, overall, the CCREM indeed provides a better representation of the correlation features than all previously proposed models.
The CCREM is constructed as a mixed model of approximately 71.5% from the ak135 model and 28.5% from the PREM2 for P-wave velocities. Similarly, it has 71.5% from the ak135 model and 28.5% from PREM in terms of S-wave velocities. Note that we approximate the velocity and density structure of ak135 and PREM using piecewise cubic polynomials in radius in this study. This results in mixed proportions in depth that are not strictly constant. Nevertheless, the CCREM is generally closer to ak135 model than to PREM. The elastic-parameter profile for CCREM is shown in Table S3. The P-wave profile of CCREM is displayed in Figure 10 as well as the S-wave profile in Figure S5 along with PREM, ak135, and ek137 models. The difference in the S-wave velocity profile in the mantle among these models is relatively small compared to the P-wave velocity profile. It is currently impossible to make more quantitative conclusions on the resolution MA AND TKALČIĆ 10.1029/2021JB022515 12 of 20 difference between the P-and S-wave velocity profiles since the exact derivation of sensitivity kernels of correlation features for Earth' velocity is not trivial.
In CCREM, we fix the crust-mantle and upper-mantle discontinuity depths (35, 410, and 660 km; corresponding to the Moho, the upper and lower discontinuities of the mantle transition zone) as in the ak135 model. Note that the differences in wave speeds among these four models are quite small except for the regions in the upper mantle and near the internal boundaries.

Travel Time Comparison
We further compare the travel times for a data set of phases through CCREM and other reference models (PREM, ak135, and ek137) in Table 2. We calculate the average travel time residual (ATTR) for each seismic phase between CCREM and another model as where N represents the number of distance data points for the specified distance range of a seismic phase with a step of 1°, and abs denotes the absolute value of travel time difference for one phase between two models. Comparisons of travel times for some seismic phases through the CCREM, PREM, ak135, and ek137 models are shown in Figure S6.
Our results show that the overall travel times through CCREM are in a much better agreement with those through short-period travel-time models (ek137 and ak135). The average value of ATTR for all phases is about 0.9 and 1.1 s for ek137 and ak135 models, respectively. However, we note that the predictions for multi-ScS phases in CCREM are closer to those in PREM. This is possibly due to the similar sensitivities of both data sets to the shear-wave structure of the mantle. Overall, travel times for the mantle phases through CCREM are closer to those through the ak135 model, while travel times for the core phases are in a better agreement with those through the ek137 model. In particular, phases with multiple P legs in the top outer core display similar ATTR values to those in the ek137 model proposed to provide a better fit to the core phases. However, most PKP-related phases show smaller ATTR values in the ak135 model than in the ek137 model.  The difference in travel times between CCREM and other models could probably arise from (a) different sensitivities for different types of data sets used (normal modes, travel times of body waves, and waveform modeling of the cross-correlation data), (b) different approaches in constructing the models and (c) different sampling coverages from data. The comparison of travel times among these models demonstrates that in general the CCREM can also provide a good representation of the travel times of regular seismic phases.

Discussion
This study proposes a new reference Earth model that provides an optimal representation of a substantial set of correlation features in global correlograms. The period band for the correlation feature data set is 15-50 s, centered between periods used in body-wave travel time and normal-mode data. Consequently, our parameterization is ineffective in resolving the thin-layered structures or relatively weak discontinuous layers inside the Earth, such as the crust, mantle transition zones, or thickness of D'' in the lowermost mantle. The sensitivity of the correlation features to these structures must be further tested within a framework of, for example, the coda-correlation wavefield formation (Wang & Tkalčić, 2020a) or full-waveform inversion (Sager et al., 2020).
Although the mantle transition-zone discontinuities in CCREM do not play a dominant role in fitting the observations, adding these discontinuities in the model slightly improves the waveform fit. Therefore, we keep the discontinuity depths as those in ak135. Initially, both the CMB and ICB in our model are fixed as those in PREM. However, the ICB depth (4 km) difference between the ak135 model and PREM is relatively significant. Compared with the CMB depth difference of only 0.5 km, this suggests that the ICB depth is much less constrained in previous reference models. We then perform synthetic tests to investigate the effects of four ICB depths (5,149.5 km (4 km shallower than ak135), 5,151.5 km (2 km shallower than ak135), 5,153.5 km (same as ak135), and 5,155.5 km (2 km deeper than ak135)) in the model on the fitness of correlation features. Our results show that models with the ICB depth of 5,155.5 km display slightly larger CC and PCC values than models with other ICB depths. However, the Wilcoxon signed-rank test shows that the mean of the CC/PCC value differences among these models is almost zero.  directly affect the travel times of the PKiKP phase, more high-quality PKiKP data are needed to strictly refine the ICB depth in the future. Here we select the ICB depth of 5,153.5 km, as in the ak135 model.
The CCREM differs only slightly from the ak135, ek137 models and PREM in the middle and lower mantle (approximately between 800 and 2,100 km) as well as the central outer core (approximately between 3,500 and 4,900 km) (Figures 10b and 10e). The significant differences among these models arise in the crust, upper mantle, and regions near the internal discontinuities, which are also the depth zones where the ak135 model differs most from PREM. The crustal model in CCREM shows a relatively larger velocity gradient compared to previous crustal models derived from travel-time data. However, the synthetic tests show that the overall fit for the whole-Earth correlation features is degraded by replacing the crustal model of CCREM with the ak135 or ek137 crust. That is probably because the two data sets (travel times of body waves vs. cross-correlation function waveforms) have different sensitivities to different structures. In the upper mantle, unlike the reduced velocity gradient between the Moho and 210 km in PREM, the CCREM shows an increased gradient in velocity and density because of the closer profile in CCREM to that in the ak135 model. However, for a spherical average, a low-velocity zone or a low-density zone is needed to match the free oscillation frequencies in PREM (Montagner & Kennett, 1996). Nonetheless, we cannot tightly constrain the crust and upper mantle structures in the current study since no specific correlation features are exclusively sensitive to these regions.
In the D'' region, the CCREM is characterized by a reduced P-wave velocity relative to PREM, ak135, and ek137 models (Figure 10c). This is because we include a portion of P-wave velocities from PREM2 D'' model in the CCREM. In PREM2, the reduced P-wave velocity in the lowermost mantle is derived as an adjustment for the separations of PKIKP and PKP ab at large distances (Song & Helmberger, 1995). Other studies using diffracted waves support a lower P-wave velocity in the D'' region relative to PREM (Garnero et al., 1993;Sylvander et al., 1997;Wysession et al., 1992). In addition, by analyzing antipodal diffracted data, Butler and Tsuboi (2020) derived a relatively lower global-mean apparent velocity within the D'' layer above the CMB. As our mid-period data suggests, if this is indeed the case, a globally averaged reduction in D'' velocities sheds important light on the chemical compositions and thermal conditions near the CMB. This aspect deserves rigorous investigation in the future.
The CCREM further shows that the P-wave velocity in the top ∼500 km of the outer core is slower, and the velocity gradient is steeper than PREM (Figure 10d). This velocity profile is more consistent with that in the outermost core of the ek137 model. Although differing in specific values of wave speed in the outer core's top, this velocity profile agrees with results from previous studies using SmKS body waves (e.g., Kaneshima & Helffrich, 2013;Tanaka, 2007;Tang et al., 2015;Wu & Irving, 2020) along with the normal mode data (van Tent et al., 2020) in a broad sense. Lower velocities than those in PREM in the outermost core (∼500 km below the CMB) could possibly imply the accumulation of light elements due to chemical reactions between the core and mantle (e.g., Buffett & Seagle, 2010) or releasing from the inner core crystallization (e.g., Franck, 1982) or primordial layering in the core (e.g., Bouffard et al., 2020). To test the validity of CCREM's outer core, comparing the predicted differential SmKS travel times with the increasing number of seismic observations is needed. However, this is beyond the scope of the current study.
In terms of the lowermost outer core, using different types of PKP data, previous studies have shown a slower P-wave profile relative to PREM (e.g., Kaneshima et al., 1994;Song & Helmberger, 1992;Souriau & Poupinet, 1991;Yu et al., 2005) and a velocity gradient generally between ak135 and PREM in the bottom ∼150 km of the outer core (e.g., Ohtaki & Kaneshima, 2015;Zou et al., 2008). It is worthwhile to note that the velocity profile in the CCREM roughly displays the same pattern. Such a P-wave velocity profile in the lowermost outer core possibly indicates a density stratification resulting from the freezing and re-melting process at the ICB (e.g., Alboussière et al., 2010;Monnereau et al., 2010;Souriau, 2015).
The Bullen parameter ŋ in the PREM approximately equals unity through the entire outer core, suggesting a homogeneous, adiabatic medium (Bullen, 1963). In comparison, the deviation of P-wave velocity in the outer core of CCREM from PREM is possibly indicative of compositional heterogeneity in the Earth's outer core (e.g., Fearn & Loper, 1981;Kaneshima & Helffrich, 2013). More detailed structures in the top and bottom parts of the outer core can be investigated by focusing on core-sensitive correlation features in future studies.
Due to strong influences from heterogeneity and anisotropy, there is a high level of uncertainty on inner core properties (Tkalčić, 2015). Moreover, most of the selected correlation features in our study are affected mainly by the mantle and outer core structures. Although some features are sensitive to the inner core structure (e.g., I2*, I4*), they are significantly affected by the mantle and outer core structures because of the small radius of the inner core. The P-wave velocity profile in the inner core of CCREM is thus loosely constrained and does not deviate much from either PREM or ak135 model (Figure 10f). In terms of the shear-wave velocities in the inner core, variable estimates have been proposed in previous studies (e.g., Cao et al., 2005;Deuss et al., 2000;Robson & Romanowicz, 2019;Tkalčić & Phạm, 2018). Based on the coda-correlation wavefield, Tkalčić and Phạm (2018) inferred an inner core model with a shear-wave velocity reduction of 2.5 ± 0.5% relative to PREM by using a single correlation feature, I2-PKJKP. This feature is more sensitive to the inner core shear-wave structure than the features selected in our study. However, we cannot casually compare the shear-wave velocities in the CCREM's inner core with results in these studies since we only focus on constraining the whole Earth's radial structure using the entire data set of correlation features. The shear-wave velocity structure in the inner core should be resolved using more specific features (e.g., I2-PKJKP, PKIKP-PKJKP) in the future.
Although seismic wave speeds play a major role in affecting the correlation features in the correlogram, we cannot avoid trade-offs between velocity and density in generating those synthesized features. Therefore, in this study, the Earth's density profile cannot be tightly constrained via fitting the correlation features. It remains ambiguous to some extent. To resolve such ambiguities in the density distribution, more normal mode eigenfrequency data can be integrated. Apart from the density distribution, Earth's radial attenuation structure will be taken into account in future coda-correlation studies, along with selected normal-mode data.
On the one hand, the inclusion of dense networks such as the USArray in this study improves the quality of global correlograms at small inter-receiver distances (the far-left side of the correlogram) because more receiver pairs exist. On the other hand, dense regional networks can introduce localized effects of the Earth's heterogeneity beneath the receivers. However, such effects will be globally averaged at larger distances since there is a benefit of the azimuthally diverse receiver-pair geometries (Tkalčić and Phạm, 2020). Also, only eight features (from the 71 used features in total) span near-zero inter-receiver distances. Therefore, the possible bias from the uneven distribution of stations on the results is minimal.
The CCREM presented here should not be considered a replacement of the ak135 or PREM models but rather a new concept based on a different data source. The ak135 model has been demonstrated to be very effective in applications for event location and predicting arrival times of various seismic phases (Kennett et al., 1995) for several decades. Meanwhile, elastic parameters in PREM have been widely used as observations compared with experimental or theoretical results in the mineral physics community. Additional observations must be considered to validate the CCREM in future studies. Hopefully, the newly proposed model could serve as a better reference model for studies using medium-period data, such as full-waveform inversion or source mechanism retrieval. Last but not least, we expect more uses of CCREM in mineral physics studies in the near future.

Conclusions
We construct the global coda correlogram by stacking cross-correlations of late-coda recordings from ten selected large earthquakes. We identify and choose a set of 71 prominent correlation features in the observed correlogram as our observations. We then derive a new spherically symmetric Earth model called CCREM, which is sensitive to the medium wave-period range. To our knowledge, this is the first reference Earth model derived from data that is not direct observations of body-wave travel times or eigenfrequencies of Earth's normal modes. Travel times and waveforms are implicit in our observables as the correlation features arise due to the similarity of seismic phases illuminating the Earth's interior after large earthquakes. The number of correlation features exceeds the seismic phases used in constructing previous Earth models because the similarity between weak seismic phases is more prominent than the weak phases themselves. CCREM is built with a combination of constraints from PREM/PREM2 and ak135 models and is designed to provide an optimal representation of the most prominent coda-correlation features in the global correlogram.
Compared with previous reference models, the CCREM displays different velocity profiles, especially in the D'' region above the CMB, the top and lowermost outer core.

Data Availability Statement
Raw seismic data are downloaded from Incorporated Research Institution for Seismology Data Management Center (IRIS DMC, https://ds.iris.edu/ds/nodes/dmc/) using SOD software (Owens et al., 2004). The information of ten events and stations used in this study is available at https://figshare.com/articles/dataset/Event-station_zip/14702745. All the figures are made with GMT6 (Wessel et al., 2019) and Matplotlib (Hunter, 2007).