Three‐Dimensional Basin and Fault Structure From a Detailed Seismic Velocity Model of Coachella Valley, Southern California

The Coachella Valley in the northern Salton Trough is known to produce destructive earthquakes, making it a high seismic hazard area. Knowledge of the seismic velocity structure and geometry of the sedimentary basins and fault zones is required to improve earthquake hazard estimates in this region. We simultaneously inverted first P wave travel times from the Southern California Seismic Network (39,998 local earthquakes) and explosions (251 land/sea shots) from the 2011 Salton Seismic Imaging Project to obtain a 3‐D seismic velocity model. Earthquakes with focal depths ≤10 km were selected to focus on the upper crustal structure. Strong lateral velocity contrasts in the top ~3 km correlate well with the surface geology, including the low‐velocity (<5 km/s) sedimentary basin and the high‐velocity crystalline basement rocks outside the valley. Sediment thickness is ~4 km in the southeastern valley near the Salton Sea and decreases to <2 km at the northwestern end of the valley. Eastward thickening of sediments toward the San Andreas fault within the valley defines Coachella Valley basin asymmetry. In the Peninsular Ranges, zones of relatively high seismic velocities (~6.4 km/s) between 2‐ and 4‐km depth may be related to Late Cretaceous mylonite rocks or older inherited basement structures. Other high‐velocity domains exist in the model down to 9‐km depth and help define crustal heterogeneity. We identify a potential fault zone in Lost Horse Valley unassociated with mapped faults in Southern California from the combined interpretation of surface geology, seismicity, and lateral velocity changes in the model.


Introduction
The Salton Trough ( Figure 1) contains active faults that make it a high earthquake hazard region in Southern California (Jennings & Bryant, 2010;Jones et al., 2008). One of the most significant threats facing this region is an overdue estimated Mw7.2 to 8.1 earthquake rupture at the southeastern end of the San Andreas fault near the Salton Sea, with energy propagating northward through the Coachella Valley into the Los Angeles basin (Day & Roten, 2012;Fialko, 2006;Jones et al., 2008;Olsen et al., 2006;Porter et al., 2011;Rose et al., 2013). Damage from such an earthquake would result in a great loss of lives and infrastructure in Southern California. To better prepare and help mitigate the risk, improved estimates of ground shaking are required. Furthermore, the importance of accurate knowledge of sedimentary basin shape and fault geometry in improving ground motion simulation estimates of seismic hazard and risk is well established (Day et al., 2008;Frankel & Vidale, 1992;Graves et al., 1998;Komatitsch et al., 2004;Passone & Mai, 2017).
In this paper, we present a 3-D seismic velocity model obtained from the simultaneous inversion of first P wave travel times from local earthquakes and explosions from the Salton Seismic Imaging Project (SSIP). Abundant seismicity along several of the main faults in the area is recorded by the Southern California Seismic Network (SCSN; Figure 2). The SSIP shots were recorded by both SSIP and SCSN stations. The combination of these two data sets provides a very dense sampling of the 3-D upper crustal seismic structure of Coachella Valley down to~9-km depth. A similar study that drew from these two data sets is Persaud et al. (2016), in which Imperial Valley to the south is modeled. Our velocity model shows strong crustal heterogeneities that correlate well with the regional geology and provides additional details on the seismic structure of fault zones. Results presented in this study should be used to improve the seismic hazard assessment for Coachella Valley by embedding our velocity model into the latest regional velocity models for Southern California.

Geological and Geophysical Setting
The Salton Trough comprises, from south to north, the Mexicali and Imperial valleys, the Salton Sea, and the Coachella Valley (Figure 1). The Salton Trough (Figure 1) is the landward continuation of a system of rightstepping transform faults in the Gulf of California separated by short incipient spreading centers (Elders et al., 1972;Larson et al., 1968;Lonsdale, 1989;Persaud et al., 2003;Stock & Hodges, 1989). However, spreading in the Salton Trough has not evolved into clear spreading centers (Elders et al., 1972;Han, Hole, Stock, Fuis, Kell, et al., 2016;Larson et al., 1968;Lonsdale, 1989;Rose et al., 2013). Sedimentary basins in the Salton Trough are located in right step overs between the Cerro Prieto and Imperial faults and the Imperial and San Andreas faults (Figure 1). Imperial Valley basin is a wedge-shaped basin, and the basin near the Cerro Prieto fault is a rhombochasm (Persaud et al., 2016). Coachella Valley contains the southern portion of the San Andreas fault system, which is the northernmost major right-lateral,strike-slip fault in the region (Lonsdale, 1989;Stock & Hodges, 1989). Geologic, gravity, and seismic studies show that Coachella Valley basin lies above a block of the Peninsular Ranges batholith that is tilted down northeastward and compressionally faulted against the San Andreas fault (Biehler, 1964;Dorsey & Langenheim, 2015;Fuis et al., 2012).
The majority of the faults in the Salton Trough are northwest striking, with similar orientations to the plate boundary ( Figure 1). However, at the northern end of the Salton Trough and to the east of Coachella Valley, there are major east-west striking faults within the North American plate including the Salt Creek fault, Smoketree Wash fault, and Blue Cut and Pinto Mountain fault zones (Figure 2). In the late Miocene to early Pleistocene, most of the displacement in this region was accommodated on the San Andreas and West Salton  Fenby and Gastil (1991), Jennings and Bryant (2010), and Rockwell et al. (2015). Pink box shows our study area in Coachella Valley (Figure 2). Black arrows in the inset map show the relative motion of the Pacific and North American plates and the red line represents the plate boundary. detachment fault systems (Dorsey & Langenheim, 2015;Jänecke et al., 2010;Mason et al., 2017;Powell, 1993;Steely et al., 2009), as such older and possibly inherited structures exist in the study area. However, in the middle Pleistocene, a series of fault systems including the San Jacinto, San Felipe, and Elsinore faults ( Figure 1) cross cut, displaced and deactivated the West Salton detachment fault (Jänecke et al., 2010;Li & Liu, 2007;Steely et al., 2009).
Orientation differences of~7°between the San Andreas fault strike and Pacific-North American slip vector results in the uplift of several hills including Durmid, Mecca, and Indio hills along the San Andreas fault system (Figure 2;Bürgmann, 1991;Fattaruso et al., 2014;Jänecke et al., 2018;Moser et al., 2017;Sylvester, 1991;Sylvester & Smith, 1976), accompanied by the formation of a flower-like structure (Fuis et al., 2017). Jänecke et al. (2018) describe a complex ladder-like structure of cross faults accommodating transpression across the San Andreas fault zone in Durmid Hill. Most of the seismicity in Coachella Valley is concentrated on the north-northwest striking faults along the northeastern side of the valley, with little seismicity located immediately west of the surface trace of the San Andreas fault (Figure 2). The rest of the seismicity in the area is clustered around the San Jacinto fault zone ( Figure 2). East-northeast striking secondary cross faults in a right step between the Coyote Creek fault and Clark fault accommodate rotating blocks in the San Jacinto fault zone . We refer to Stock and Hodges (1989) for more information on the tectonic evolution of the Salton Trough region prior to 6 Ma.  Hauksson et al. (2012) catalog, selected to a maximum depth of 10 km. Light blue triangles are the SCSN stations. Black lines are surface traces of mapped faults in the area from Jennings and Bryant (2010). Yellow lines are gravity-derived basement contours labeled every 1 km from Dorsey and Langenheim (2015) and Langenheim et al. (2005) The Coachella Valley basin is filled with late Cenozoic sedimentary rocks, deposited from the adjacent mountain ranges and the Colorado River Dorsey, 2012;McNabb et al., 2017;Robinson et al., 1976;Sylvester & Smith, 1976). The northwest trending Peninsular Ranges batholith and Eastern Transverse Ranges located on the western and eastern sides of the trough, respectively, are mostly composed of Mesozoic crystalline basement rocks McNabb et al., 2017;Morton et al., 2014;Powell, 1993;Sylvester & Smith, 1976). Rocks in the Peninsular Ranges are mostly plutonic, comprising granodiorite and tonalite, compared to the Eastern Transverse Ranges rocks, which are more diverse and include volcanic and metamorphic rocks ( Figure 3). The Peninsular Ranges batholith is compositionally divided into eastern and western provinces Magistrale & Sanders, 1995 and references therein;Morton et al., 2014). The Western Peninsular Ranges batholith is older than the eastern section, and unlike the plutons in the western section that include gabbro, the eastern section has little to no gabbro (Magistrale & Sanders, 1995, and references therein). There are zones of mylonitic rocks resulting from highly deformed granodioritic and tonalitic protoliths that have been mapped in the Eastern Peninsular Ranges, adjacent to the Salton Trough, known as the Eastern Peninsular Ranges mylonite zone (EPRMZ). Kern and Wenk (1990) performed seismic laboratory experiments on rocks from the EPRMZ at different temperatures and pressures similar to conditions at which the rocks may be found at depth and report velocity anisotropy of~1.7% to 19%. Allam and Ben-Zion (2012) and Fang et al. (2016) provide evidence from their seismic velocity models for displacement of the Peninsular Ranges across the San Jacinto fault zone. Sylvester and Smith (1976) and McNabb et al. (2017) describe, in detail, the lithology and structural features of exposed, pre-Cenozoic basement rocks between the Skeleton Canyon and Painted Canyon faults in the Mecca Hills.

Previous Geophysical Studies
Several regional tomographic models have been developed for Southern California and incorporated into physics-based hazard analysis tools to forecast ground shaking. Chen and Lee (2015) summarize the historical development of these models. Southern California Earthquake Center (SCEC) researchers have developed two regional velocity models, including the Community Velocity Model-Harvard (CVM-H) and Community Velocity Model-SCEC (CVM-S). Previous versions of CVM-S and CVM-H included a regional background seismic velocity model from Hauksson (2000) and Magistrale et al. (2000) travel time tomography. CVM-H includes a near-surface geotechnical layer from Ely et al. (2010), crust-mantle interface from Yan and Clayton (2007) receiver function studies, upper mantle velocity model from Prindle and Tanimoto (2006) finite-frequency tomography, and high-resolution seismic velocity models of Southern California basins, developed by Süss and Shaw (2003), using oil industry well logs and seismic reflection data. CVM-H 15.1.0 is the latest version of CVM-H and is similar to CVM-H 11.9, created by improving the previous versions, using full 3-D waveform inversion (Tape et al., 2009(Tape et al., , 2010. Barak et al. (2015) used ambient noise seismic tomography to improve CVM-H 11.9. CVM-S also includes a geotechnical layer, crust-mantle interface from Zhu and Kanamori (2000), upper-mantle velocity structure from Kohler et al. (2003), and embedded, detailed basin velocity models from Magistrale et al. (1996), which used a velocitydepth-age function determined from geological data. The latest version of CVM-S is CVM-S 4.26, which was developed using the previous version (CVM-S 4) as a starting velocity model in a full 3-D tomographic inversion (Lee, Chen, Jordan, Maechling, et al., 2014). The latest CVM-H and CVM-S models have been shown to improve calculated waveform misfits with observed seismograms, compared to their predecessors (Chen & Lee, 2015;. Different 2-D and 3-D local seismic velocity models have been developed for several parts of the Salton Trough. However, until recently, results from these studies focused more on Imperial Valley to the south and have limited discussions of Coachella Valley, as they lack dense data coverage (Biehler, 1964;Fuis et al., 1982Fuis et al., , 1984Parsons & McCarthy, 1996). The SSIP (Figure 1), which consists of refraction and wideangle reflection profiles, was carried out in early 2011 to primarily investigate earthquake hazard in the Salton Trough region (Rose et al., 2013). SSIP results for 2-D profiles across the Coachella and Imperial valleys, Salton Sea, and along the axis of the Salton Trough have been analyzed (Fuis et al., 2017;Han, Hole, Stock, Fuis, Kell, et al., 2016;Han, Hole, Stock, Fuis, Williams, et al., 2016;Hernandez et al., 2015;Persaud et al., 2014;Persaud et al., 2016;Sahakian et al., 2016). Fuis et al. (2017) focused on imaging the subsurface geometry of the San Andreas fault system and other faults in Coachella Valley along SSIP Lines 4 and 6 ( Figure 2), by migrating reverse-moveout phases using the method of Bauer et al. (2013). Nicholson (1996) had previously analyzed aftershock sequences and focal mechanisms from the 1948 Desert Hot Springs and 1986 North Palm Springs earthquakes to infer a curvilinear geometry for the Banning strand of the San Andreas fault. Fuis et al. (2012) develop a model of a nonvertical San Andreas fault through much of Southern California, including the Coachella Valley section where the fault dips to the northeast, using seismicity, potential field data, and seismic tomography interpretations. Han, Hole, Stock, Fuis, Kell, et al., (2016) produced a seismic velocity model for SSIP Line 1 (Figure 1), along the axis of the Salton Trough. Dorsey and Langenheim (2015) interpreted the isostatic gravity map of Coachella Valley and analyzed a gravity transect across the valley just north of the Salton Sea close to SSIP Line 4. Lin (2013) inverted P wave and S wave travel times from more than 132,000 local earthquakes using the simul2000 tomography algorithm (Eberhart-Phillips, 1990;Thurber, 1983Thurber, , 1993Thurber & Eberhart-Phillips,1999) to produce a 3-DVp and Vp/Vs velocity model for the Salton Trough. Their model has a horizontal grid spacing of 5 km and a variable vertical grid spacing between 2 and 4 km. Final hypocenter locations from the Lin (2013) inversion yielded a tighter seismicity pattern than the SCSN catalog in the area that correlates well with mapped faults.

Travel Time Data
Our study area and the seismic array that recorded the data used in our tomographic inversion are shown in Figure 2. We utilized only first arrival times, which include mostly Pg and direct or shallowly refracted phases because of our interest in the upper crustal structure. Prior SSIP studies such as Han, Hole, Stock, Fuis, Kell, et al. (2016), Han, Hole, Stock, Fuis, Williams, et al. (2016) used later arrivals to determine deeper crustal structures. There has also been work that invert amplitude data from regional seismic networks to obtain an attenuation model for the Salton Trough, for example, Lin (2014); other studies like Lee, Chen, Jordan, Maechling, et al. (2014), Tape et al. (2009), andTape et al. (2010) utilize entire seismograms in a full 3-D waveform inversion to develop regional velocity models of Southern California. 3.1.1. Salton Seismic Imaging Project Data A major part of our data set is from SSIP, which consists of a series of intersecting refraction and wide-angle reflection profiles, and 2-D receiver arrays with dense grid spacing across the Salton Trough (Han, Hole, Stock, Fuis, Kell, et al., 2016;Rose et al., 2013). Two types of portable land seismometers were deployed in the SSIP: Ref Tek 130's ("RT130"), with 4.5-Hz three-component sensors and Ref Tek 125a's (Texans), with 4.5-Hz single-component sensors (Rose et al., 2013). Lines 4 and 7 also had cabled geophones at 10-m spacing for a couple of kilometers, spanning the San Andreas fault. Ocean bottom seismometers deployed in the Salton Sea were L-CHEAPO 4 × 4's, with 4.5-Hz four-component sensors (Rose et al., 2013). Land shots were carefully designed to prevent damage in surrounding areas and were made by detonating nitrate-based explosive agents in deep boreholes, where the center of the charge ranged from about 20-to 30-m depth, with an equivalent energy release of magnitude 1 to 2 earthquakes (Rose et al., 2013). Air gun shots were fired at 2-to 4-m depth in the Salton Sea. In summary, the SSIP comprises 2,709 seismometers deployed at 4,341 sites on land, 48 ocean bottom seismometers deployed at 78 sites in the Salton Sea, 126 land shots, and 2381 air gun shots (Rose et al., 2013). Refraction data from the SSIP survey used in the inversion include 44,823 first P wave travel times that were manually picked from 251 land and sea shots recorded at 1,919 SSIP stations and 82 SCSN stations. Shots outside Coachella Valley ( Figure 1) were included wherever possible. The uncertainty in the picks depends on the shot size and proximity to receivers (Figures 1 and 2). SSIP shots that are farther away from the areal extent in Figure 1 (pink box in Figure 2) or shots with smaller explosives weight tend to be associated with larger root-mean-square (RMS) travel time residuals ( Figure S1 in the supporting information). Han, Hole, Stock, Fuis, Kell, et al. (2016) and Han, Hole, Stock, Fuis, Williams, et al. (2016) give discussions on the SSIP data and provide estimates of pick uncertainty.

Hauksson et al. (2012) Hypocenters
In addition, earthquake travel time data from events in the Hauksson et al. (2012) catalog were obtained from the Southern California Earthquake Data Center using the Seismogram Transfer Program (Southern California Earthquake Data Center, 2013; http://scedc.caltech.edu/research-tools/stp-index.html). The earthquake catalog comprises events recorded by the SCSN from 1981 to June 2011 and relocated from their initial hypocenters, determined using a 1-D velocity model. Hauksson et al. (2012) relocate the hypocenters in a two-step procedure. A 3-D velocity model from Hauksson (2000) is first used to obtain more accurate absolute earthquake locations than the initial catalog hypocenters. Hypocenters relocated using this 3-D velocity model are more accurate than a 1-D velocity model, which lacks lateral heterogeneity (Thurber, 1983). The second step utilizes waveform cross correlations Lin et al., 2007;Shearer et al., 2005) between similar events to improve the relative hypocenter locations. The relative earthquake locations usually exhibit a tighter spatial clustering that aligns well with Quaternary faults and is 10.1029/2018JB016260 Journal of Geophysical Research: Solid Earth more relatable to the regional geology than absolute earthquake locations Lin et al., 2007).
Local events in the catalog were selected to a maximum focal depth of 10 km to focus on the upper crustal structure. Travel times for events in the catalog were not included in the inversion if the event was not recorded by at least four stations located in our study area. In the earthquake data set used in our study, 88% of the events were relocated using cross correlation, and~11.9% of the earthquake locations were determined using the 3-D velocity model from Hauksson (2000). Figure 2 shows the 39,998 earthquakes recorded by 76 SCSN stations used in the inversion, resulting in a total of 372,236 earthquake travel times. Most of the earthquakes are concentrated on the San Jacinto fault zone in the Peninsular Ranges and on faults at the northwestern end of the Little San Bernardino Mountains in Figure 2. The spatial resolution of the inversion is controlled mainly by the angular coverage of the crossing source-receiver paths. Shots on each SSIP line were recorded by stations on other SSIP lines and on the areal grid to provide a dense crustal sampling (Rose et al., 2013). The SSIP data complement the earthquakes by providing coverage within the valley and of the shallow velocity structure, which are usually lacking in earthquake-only tomographic studies.

Tomographic Inversion Method
Our preferred starting velocity model is the 1-D average of the 3-D velocity model from Lin (2013) in Figure 4a. We also tested the 1-D averaged Hadley and Kanamori (1977) and the SCEC CVM-H 11.9.1 models shown in Figure 4a, which produced similar features in the final velocity model. We modified the starting velocity models by extending them up to 5-km elevation using linear velocity-depth relationships.

Simultaneous Inversion Technique
A finite difference scheme for the eikonal equation (Hole & Zelt, 1995;Vidale, 1990) is first used to compute the travel times in the starting velocity model. We then iteratively invert the travel time data for velocity model perturbations using the nonlinear travel time tomography algorithm from Hole (1992) and Hole et al. (2000), which allows for hypocenter inversion. The Hole et al. (2000) hypocenter relocation algorithm uses singular value decomposition to invert for changes in the hypocenter locations that lead to a reduction in the observed travel time misfits. Equations of the objective functions for hypocenter relocations in the Hole et al. (2000) algorithm and SIMULPS (Eberhart-Phillips, 1990;Thurber, 1992Thurber, , 1993 or HYPOINVERSE (Klein, 2002) programs used in Hauksson et al. (2012) are similar. However, the latter programs use a different technique to invert for the location perturbation matrices. Also, unlike the SIMULPS and HYPOINVERSE programs that calculate travel times using ray tracing, the Hole et al. (2000) algorithm uses a finite difference approach (Vidale, 1990). As such, the Hole et al. (2000) algorithm, like SIMULPS or HYPOINVERSE, is only capable of improving absolute earthquake locations, especially when the original catalog hypocenters are determined using a 1-D velocity model. Furthermore, this approach does not guarantee a tighter clustering of hypocenters, because the hypocenters move independently from each other and can become more spread out from their initial or actual locations rather than focused.
The inversion strategy is based on the minimum structure approach (Hole et al., 2000;Hole et al., 2006). We experimented with different grid spacings and gridded our velocity model at 0.5-km fixed spacing within a 214 by 198-km area to ensure a suitable level of detail based on the distribution and density of crossing raypaths. Our model area is larger than what is shown in Figure 2 because we include SSIP shots in Imperial Valley. Depths are given with negative values above sea level and range from −5 to 16 km. Topography is included in the velocity model so that station time corrections are not used. The smoothing dimensions of the model perturbations were changed after every three tomographic iterations. Large smoothing operators of 300 by 300 by 30 grid cells were used in the early iteration to model the large-scale features. The smoothing size was reduced to 10 by 10 by 2 grid cells in the later iterations to reveal small-scale structures, thus ensuring a slow, stable, and global minimum convergence (Figure 4b).

Preferred Inversion Scheme
In local earthquake tomography, it is generally considered incorrect to invert for model parameters without hypocenter relocations (Hole et al., 2000;Thurber, 1992), because the original catalog hypocenter locations are deemed inaccurate. In our study, we simultaneously invert explosion and earthquake travel times and initially experimented with relocating the Hauksson et al. (2012) hypocenters, but as explained below, our chosen inversion scheme uses fixed hypocenters.

10.1029/2018JB016260
Journal of Geophysical Research: Solid Earth  (2013) is our preferred model. We modified the starting velocity models by extending them up to 5-km elevation using linear velocity-depth relationships. (b) RMS travel time residuals for each iteration in the tomographic inversions using three starting velocity models. The residuals decrease in a similar and stable manner for all three starting models. Gray star represents the RMS residual for the first velocity inversion iteration. The 26th iteration marked with blue star in the Lin (

Journal of Geophysical Research: Solid Earth
We tested different inversion schemes using different numbers of interleaving tomographic and hypocenter iterations and varying weights of the SSIP explosion data (Table S1). A similar smoothing scheme is employed for each set of tomographic and hypocenter iterations to ensure that the model converged at the same tomographic iteration number for all of our experiments, which also facilitated the comparison of the final model results. In general, the relocated hypocenter locations from our different inversion schemes were similar to the Hauksson et al. (2012) hypocenters. To summarize the results presented in Table S1, the earlier the hypocenters were relocated in the inversion, the less clustered the results were relative to the later relocated hypocenters. Moreover, the more hypocenter relocations performed after the tomographic iterations, led to less clustered hypocenter locations. None of our inversion schemes that inverted for the hypocenter locations produced more tightly clustered hypocenters than the Hauksson et al. (2012) catalog. An important conclusion from these tests is in the case of our study; it is better to fix the hypocenter locations throughout the inversion rather than relocate a catalog that was produced through an extra step of waveform cross correlation, which might introduce errors into our model. Furthermore, because the earthquake catalog relates well with the geology, we have high confidence in regions of our model constrained only by the hypocenters. Our final velocity model, chosen from the 26th velocity iteration (blue star in Figure 4b), has an RMS residual of 199 ms and is free from modeling artifacts such as the bull's-eye artifact resulting from noise overfitting. There was also a systematic and stable reduction in the travel time residuals in the tomographic iterations leading up to the 26th, and subsequent iterations had similar residual distributions ( Figure S2). RMS residuals from the 26th iteration for the inversions that use a 1-D-averaged Hadley and Kanamori (1977) model and SCEC CVM-H 11.9.1 are 214 and 191 ms, respectively. Histogram distributions comparing the 1st and 26th iteration for the preferred inversion show the model improvement, with the final travel time residuals having a strong peak centered about zero (Figure 4c).

Model Resolution Evaluation
We employ a combination of ray coverage maps (Figures 5 and S3) and synthetic model reconstruction tests (Figures 6 and S4-S8) to assess the resolution of our final velocity model. Ray coverage maps are shown for 2-,3-,6-, and 9-km depths in Figure 5. We show ray coverage for other depths in Figure S3. Above 2-km depth, ray coverage is limited in areas away from the SSIP and SCSN receivers because distant rays become near vertical closer to the stations and only sample regions directly beneath the receivers in the shallow layers. Ray coverage deteriorates rapidly below 9-km depth due to the lack of earthquakes based on our selection criterion and limited long-offset SSIP picks. Synthetic model reconstructions were used to qualitatively estimate the regional spatial resolution of our velocity model. The synthetic models, similar to checkerboard models, were created by adding sine function perturbations of varying wavelengths to a modified Lin (2013) 1-D-averaged starting velocity model to produce a checkerboard pattern in the velocities, with a maximum P wave velocity perturbation of 5%. Travel times were computed and inverted through the synthetic models using the same data and inversion parameters and the same technique as our actual inversion in section 3.2.1. Recovered structures show that our model can resolve features as small as 6 by 6 by 1.5 km. However, this does not imply that we cannot resolve smaller structures, especially along the SSIP Lines where ray coverage is dense. Larger structures are resolved better down to 9-km depth in areas with good ray coverage. We recover northwest trending anomalies ( Figure S7), with similar strikes as the mountain ranges and Coachella Valley to qualitatively assess the preservation of linear features in our model (Lévěque et al., 1993). For example, this allows us to see how much a continuous valley would be broken up into subbasins. In addition, seismic velocities in the selected cross sections from our velocity model are consistent with results from previous SSIP studies along Lines 1N (Han, Hole, Stock, Fuis, Kell, et al., 2016), 4 and 6 (Fuis et al., 2017), and 5 , which increases our confidence in other areas of our 3-D velocity model.

Results
Together with the ray coverage maps (Figures 5 and S3), the accompanying synthetic resolution tests (Figures 6 and S4-S8) provide model assessment information, by illuminating regions and features in our velocity model that are reliable enough for discussion. We present depth slices and vertical profiles of our 10.1029/2018JB016260 Journal of Geophysical Research: Solid Earth model in Figures 7-9 and show comparisons to previous models in Figure 10; we do not show velocities in areas with no ray coverage.
The spatial extent of good ray coverage is restricted at 2-and 3-km depth (Figures 5 and S3). Where ray coverage is good, the average velocity in the adjacent mountain ranges is~6 km/s, while velocities in the valley range from 3 to 5 km/s (Figure 7). Synthetic resolution tests at these depths (Figures 6 and S4-8) recover the checkerboard pattern, with similar amplitudes as the input model velocities within areas with good ray coverage. However, due to limited ray coverage away from the SSIP and SCSN receivers at 2-km depth, the southeastern end of the Santa Rosa Mountains, northwestern end of the San Jacinto Mountains, and the western part of the Transverse Ranges are not imaged (white areas in Figure 7a). The coverage improves at 3-km depth and below, especially at the northern end of the San Jacinto Mountains. However, the western side of the Transverse Ranges and the southeastern end of the Santa Rosa Mountains are still poorly imaged at 3-km depth (Figure 7b).
Ray coverage is best at~4-to 7-km depths (Figures 5 and S3), where most of the valley and mountain ranges are imaged (Figures 7 and 8). These depths have the best coverage because of the intersection of turning rays from the SSIP data and upgoing refracted rays from the earthquakes, which sample a broad region in our Figure 5. Ray coverage maps for our preferred model iteration at four different depths, 2, 3, 6, and 9 km. We show ray coverage for other depths in Figure S3. White and gray areas are regions illuminated by rays. The edges of our model shown in light brown have lower coverage due to a lack of earthquakes and sparse seismic stations in those areas. Coverage deteriorates rapidly below 9 km, where we limit our discussion of the velocity model. Blue circles are Salton Seismic Imaging Project receivers and orange triangles are Southern California Seismic Network stations.

10.1029/2018JB016260
Journal of Geophysical Research: Solid Earth model. Resolution tests at these depths (e.g.,6-km depth in Figures 6 and S4-S8) show good recovery of the velocity structure of the input models in areas with good ray coverage.
Below 7-km depth, ray coverage remains good where there are a lot of earthquakes, for example, at the northeastern side of Coachella Valley and in the San Jacinto fault zone (Figure 8). The velocity structure is well imaged in the synthetic resolution tests. However, ray smearing in the horizontal direction along dominant raypaths with limited azimuths (e.g., raypaths at 9-km depth in Figure 5) becomes more pronounced at these depths. Below 9-km depth, ray coverage is strongly reduced due to our earthquake depth selection criterion, so our velocity model is not presented. Figure 10 shows cross-sectional profiles of our P wave velocities compared with CVM-H 15.1.0 (Tape et al., 2009(Tape et al., , 2010 and CVM-S 4.26 (Lee, Chen, Jordan, Maechling, et al., 2014). We assess the model differences by giving a qualitative description of the models in terms of the basin geometry and structure, basement topographic structure, and crustal heterogeneities. There are significant differences between the models, which could result in substantial changes in ground shaking hazard estimates in the region. For example, the symmetric basin in the CVMs versus the asymmetric basin defined in our model or the absence of strong upper crustal heterogeneity in CVM-H 15.1.0. To fully evaluate the model differences described in this section, we would need to embed our upper crustal P wave velocity model in either of the latest Figure 6. Synthetic resolution test for our preferred inversion parameters at 3-, 6-, and 9-km depths using 10 by 10 by 3 km checkers, with a maximum P wave velocity perturbation of 5% from Lin (2013) 1-D-averaged starting velocity model. Black lines on the 6-km depth slice in the top panel show the locations of the profiles below. Labels at the bottom right of each plot indicate input versus recovered models. In the vertical sections, information about the orientation and vertical exaggeration is plotted at the top and bottom, respectively. Black circles are Salton Seismic Imaging Project receivers, yellow squares are SSIP shots, and orange triangles are Southern California Seismic Network stations.

Journal of Geophysical Research: Solid Earth
CVMs and simulate wave propagations to compare synthetic seismograms with observed recordings of the simulated earthquakes.

Basin and Basement Topography
In the profiles shown, we highlight the shape of the basin and basement topography using the 3-to 5-km/sP wave velocity contours (Figure 10). Based on velocities of~4 km/s recorded at shot points located on bedrock  Fuis et al. (2017) report that the minimum depth to basement in the area is probably at or below the 4-km/s velocity contour. Also, Magistrale and Sanders (1995) note that the pervasive cracking of rocks might reduce the observed velocities at the surface due to a change in the mechanical properties of the rocks. For this reason, we infer that the basement surface in the area probably lies between the 4-and 5-km/s The basin geometry and basement topography in all three models are different. Along SSIP Line 1N, the sedimentary basin thickness in the models is greatest at the southeastern end and thins to the north (Figures 9a  and 10). However, in CVM-H 15.1.0, the maximum sedimentary thickness is less than in our model by a few hundreds of meters, and the sedimentary thickness reduces much quicker to the northwest, about halfway through the profile compared to our model where the sedimentary layer thickness is~2 km at the northwestern end. The basement topography represented by CVM-H 15.1.0 is more uniform and increasingly shallower to the northwest, whereas undulations in the basement surface in our model suggest the presence of multiple subbasins. CVM-S 4.26 has a much closer resemblance to our model in terms of the basin geometry and basement topography. For example, at~30 km in CVM-S 4.26, there appears to be a basement high similar to what we observe around~35-km profile distance in our model, which supports the presence of subbasins. However, seismic velocities within the basin in CVM-S 4.26 are much lower (<3 km/s) than in our model and the transition to the basement is much more abrupt.
Along SSIP Line 4, the basin geometry represented in CVM-H 15.1.0 is markedly different from our model. The maximum sedimentary layer thickness is lower, and the thickness is uniform in the southwestern end before starting to thin about halfway in the profile toward the northeast, whereas our basin geometry appears asymmetric. Again, the basin geometry in CVM-S 4.26 closely resembles our model, but the basin appears to be more symmetric. Velocities within the basin are also lower than in our model.
In SSIP Line 5, the basin geometry and basement topography in CVM-H 15.1.0 and CVM-S 4.26 are identical. The basin geometries are symmetrical in comparison to our model, which is asymmetric. Also, in both Figure 9. Cross-sectional slices through our velocity model oriented along Salton Seismic Imaging Project Lines 1N, 4, 5, and 6. Contour intervals of 1 km/s are used to highlight the shape of the sedimentary basin. Velocities in areas with no ray coverage are not shown. Yellow dots represent projected earthquakes within ±10 km from our profiles, which do not always occur on faults that cross the profiles, for example, at the northeast end of Line 6. The comparison of our model with the latest Community Velocity Models is shown in Figure 10. Information about the orientation and vertical exaggeration of the profiles is plotted at the top and bottom, respectively. White stars are projected shots from SSIP. CVH = Coachella Valley high; PRH = Peninsular Ranges high; OH = Orocopia high. Along SSIP Line 6, the sedimentary basin is thin and almost absent in the CVMs. CVM-S 4.26 has a thicker and wider sedimentary basin that is more comparable to our model. Velocities in the basin are also much lower in the CVMs than in our model.

Crustal Heterogeneities
We use the 6.4-km/s velocity contour to indicate regions of higher crustal velocities in the model than an average 6-km/s crust ( Figure 10). Along SSIP Lines 1N, 4, 5, and 6 in CVM-H 15.1.0, the crustal velocities are uniform at~6 km/s, and strong upper crustal heterogeneities are absent compared to CVM-S 4.26 and our model where there are regions with velocities that exceed 6.4 km/s. CVM-S 4.26 shows similar trends in crustal heterogeneities compared to our model. However, these high velocity areas are more discontinuous and depth restricted in our model.
Crustal seismic velocities in CVM-S 4.26 increase to >6.4 km/s around 5-to 6-km depth in all profiles in Figure 10. In SSIP Line 1N, the trend of the 6.4-km/s velocity contour in CVM-S 4.26 is similar to the top of the high-velocity areas in our model. At the southeastern end of SSIP Line 1N in our model and CVM-S 4.26, the top of this zone dips downward to the northwest to~8-km depth at~20 km along the profiles, then rises to~5-km depth at~30 km distance and remains relatively horizontal to the northwest. The highvelocity region in SSIP Line 4 is uniform in CVM-S 4.26 compared to our model where it is more discontinuous. Along SSIP Line 5 in CVM-S 4.26, the 6.4-km/s contour is dipping downward in the northeastern direction similar to the high-velocity areas that we observe in our model. At the southwestern end of SSIP Line 6, we observe high crustal velocities in both our model and CVM-S 4.26. However, the high-velocity region is much broader and extends to more than half of the profile's length toward the northeast in CVM-S 4.26.

Discussion
We discuss our velocity model shown in Figures 7-9 by associating high and low velocities, and lateral velocity changes in the model with the surface geology (Figure 3), surface traces of mapped faults in the area, and seismicity trends. Indeed, without independent constraints from other data such as well data, there is no absolute certainty in the correlations.

Structural Features in the Upper Crust
Lateral seismic velocity contrasts are strongest at shallow depths (down to~3 km), where velocities have the largest range (Figures 7a and 7b). The near-surface (2 km in Figure 7a) velocity structure correlates well with the surface geology ( Figure 3) and the isostatic gravity field in Dorsey and Langenheim (2015) and Langenheim et al. (2005). We focus our discussion of the near-surface starting at 2-km depth due to limited ray coverage at shallower depths, away from the SSIP array and SCSN stations. In general, the 2-and 3-km depth slice maps are very similar. Velocities of sedimentary rocks within the valley range from~1.5 km/s to just over 4 km/s, with a similar geometry to the Quaternary alluvium mapped at the surface ( Figure 3). Lower velocities are observed on the eastern side of Coachella Valley basin close to the San Andreas fault, hinting at the basin asymmetry (e.g., in Figure 7b). The sediments in the valley appear to be separated into subbasins by a major northwest step-up in the basement (Figure 7b), associated with a northeast trending ridge near Point Happy (PH, Figures 3 and 7b) that extends into the valley, across from the junction of the Banning and Mission creek strands of the San Andreas fault (Dorsey & Langenheim, 2015). Lower velocities (<5 km/s) persist to~4-km depth near the northern end of the Salton Sea (Figure 7d).
The Peninsular Ranges are associated with velocities that locally exceed 6.4 km/s. These zones of relatively high velocities (>6.4 km/s) within an average background velocity of~6 km/s in the Peninsular Ranges exist between 2-and 4-km depth (e.g., Peninsular Ranges high (PRH) in Figures 7 and 9). PRH could be possibly related to the anisotropic mylonitic rocks of the Santa Rosa Mountains (Erskine & Wenk, 1985;Kern & Wenk, 1990;Sharp, 1979). Our inversion method makes no assumption about anisotropy; hence, velocities reported are average velocities. However, azimuthal velocity variations of mylonitic rocks ranging from 1.7% up to 19% can be expected based on laboratory measurements by Kern and Wenk (1990). Arrival time versus azimuth analysis can be performed, and P wave azimuthal slowness perturbations can be incorporated into the inversion to study the Eastern Peninsular Ranges mylonite zone anisotropy, in detail (Kuo- Chen et al., 2013;Środa, 2006). Velocities beneath the Orocopia Mountains are poorly determined at 2-and 3-km depths, and the southeastern end of the Santa Rosa Mountains near the western edge of the Salton Sea also lacks sufficient ray coverage due to relatively sparse seismic stations in those areas (Figures 7a and 7b).

Journal of Geophysical Research: Solid Earth
Between 3-and 5-km depths, we observe low velocities beneath the northwest end of Coachella Valley basin by the Indio Hills (Coachella Valley low -CVL, Figures 7b-7d). The Coachella Valley low appears to be bounded to the southwest and northeast by the Banning and Mission Creek faults (Nicholson & Lees, 1992). Comparing the locations of the Coachella Valley Low at 3-and 5-km depths implies that the Banning and Mission Creek faults are northeast dipping (Nicholson, 1996;Nicholson & Lees, 1992) and are not a single fault plane at ≤5-km depth in this region. We note, however, that this interpretation does not suggest that the faults necessarily merge into a single fault plane at >5-km depth.
At 4-and 5-km depths (Figures 7c and 7d), we observe high velocities, exceeding 6.4 km/s, beneath the Orocopia mountains (Orocopia high, and at the southeast end of the Santa Rosa Mountains. Here the Orocopia high might be related to schistose rocks in the Orocopia mountains ( Figure 3). At~9-km depth, the maximum depth constrained by our model, the high velocities of the Orocopia Mountains extend northwest of the northeast end of SSIP Line 4 relative to their previous location to the southeast in the 5-km slice (e.g., in Figures 7c, 7d, and 8), suggesting that OH probably deepens toward the northwest. These high velocities are bounded to the north by the Smoketree Wash fault (in Figure 8).
Below~5-km depth, high velocities of~6.3 to 6.4 km/s (Coachella Valley high,CVH, in the Peninsular Ranges batholith exist beneath the thickest part of Coachella Valley basin and can be clearly seen in SSIP Lines 1N, 4, and 5 ( Figure 9). The Eastern Peninsular Ranges mylonite zone is reported to be dipping downward to the northeast (e.g., Erskine & Wenk, 1985), for this reason, we hypothesize that CVH is likely related to PRH (e.g., Figure 9c) and could be related to mylonitic rocks at greater depths beneath the valley. We attribute the lower velocities of CVH (Figures 8 and 9) beneath the valley compared to PRH observed at shallower depths to smearing along raypaths, which implies that CVH may be more depth restricted (<2 km in thickness) than it appears in the cross-sectional velocity model (Figure 9).
At 8-km depth, we infer the subsurface location of the San Andreas fault by projecting the surface trace using a 50°northeast dip (white dashed line in Figure 8c). Here the San Andreas fault appears to separate Coachella Valley high from Orocopia high (e.g., Figure 9b), which are features of the Peninsular Ranges and Eastern Transverse Ranges, respectively. The complex relationship expressed at this depth might suggest that the San Andreas fault could be coincident with the boundary of two different basement types in the region. Although we have correlated high-velocity regions (CVH, OH, and PRH) in our velocity model to geologic units mapped at the surface for discussion purposes, some of these features could also represent older inherited structures. Nonetheless, the strong crustal heterogeneities present in our model can help with modeling the rheological properties of the crust and their interaction with main faults in the area.

San Andreas and San Jacinto Fault Zones Structure
Abrupt lateral velocity changes exist across the major faults in the area, especially in the San Jacinto and San Andreas fault zones (Figures 7 and 8). We observe these lateral velocity changes in the region below the surface trace of the faults to 9-km depth. Between 2-and 5-km depths, the San Andreas fault is located on the eastern flank of the Coachella Valley basin and separates lower velocity sedimentary rocks from highervelocity sedimentary rocks (Figure 7). At greater depths of ≥6 km, the San Andreas fault likely coincides with the boundary between the Peninsular Ranges block and the Eastern Transverse Ranges (white line in Figure 8c). On the western side of the valley, at 2-to 3-km depth, velocity changes exist across the San Jacinto fault separating lower velocity (~5.6 km/s) Western Peninsular Ranges rocks from higher velocity (6 to 6.4 km/s) Eastern Peninsular Ranges rocks (Figures 7a and 7b). At greater depths, this change disappears and the SJFZ appears to be associated chiefly with high-velocity bodies.
Lateral velocity changes exist across the Smoketree Wash fault where, at and below 5-km depth, it appears to separate high-velocity rocks to the south from the~6-km/s basement rocks to the north (Figures 7 and 8). In the 7-km depth slice, we observe relatively low-velocity zones (~5.9 km/s) outlined by the 6-km/s contours (Figure 8b) that could represent apparent offset across the left-lateral Blue Cut fault. This apparent offset is roughly 10 km, measured from the centers of relatively low-velocity zones, and is less evident at shallower depths. Prior estimates of offset across the Blue Cut fault zone from geologic data (Hope, 1966) and magnetic anomalies (Langenheim & Powell, 2009) range from 3 to 9 km. Lateral velocity changes also exist from north to south across the Pinto Mountain fault zone and become more evident at ≥4-km depth (Figures 7c, 7d, 8a, and 8b).

10.1029/2018JB016260
Journal of Geophysical Research: Solid Earth

Seismicity Trends in the Mountain Ranges
Recent earthquakes in the study area (Figures 2, 7, and 8) align reasonably well with surface traces of mapped faults in the area and are mainly concentrated along the north striking faults between the Blue Cut and Pinto Mountain fault zones and along fault branches within the San Jacinto fault zone. Black dots in Figures 7 and 8 show projected hypocenter locations from the Hauksson et al. (2012) catalog located within ±500 m from the depth slices. Most of the seismicity starts at ≥2-km depth, with the majority of the earthquakes located at ≥6-km depth. At ≥4-km depth, we observe a northeast trending pattern in the seismicity that initiates at the southern end of the north-south trending faults, south of the Pinto Mountain fault zone, and extends northeastward toward the Pinto Mountain fault zone (e.g., in Figures 2  and 8b). This northeast trend had previously been identified by Hauksson et al. (1993) and is mostly related to foreshocks and aftershocks from the 1992 Joshua Tree earthquake ( Figure S9). Focal mechanism solutions show that majority of these earthquakes had a right lateral sense of slip (Hauksson et al., 1993;Yang et al., 2012). In addition to the seismicity trend, we also observe lateral velocity changes across this zone at 5-to 7km depths (Figures 7d, 8a, and 8b). However, these velocity changes are less apparent at shallower (≤4 km) and deeper (≥8 km) depths. The surface location of the velocity change and seismicity trend corresponds to Lost Horse Valley, which can be seen in the more detailed geologic map of the area (Dibblee, 1968; https:// ngmdb.usgs.gov/Prodesc/proddesc_112.htm). The valley can also be identified by looking closely at our topographic map in Figure 2. We estimate~3-km right-lateral offset based on observed displacement of gneissic rock units across the Lost Horse Valley in the geologic map of Dibblee (1968). Based on the seismicity trend, lateral velocity change in our model across this seismicity zone, and our geologic observations of displacement, it is likely the hypocenters are related to an unmapped fault in the Little San Bernardino Mountains that we name the Lost Horse Valley fault zone, shown with a white dashed line in Figure 8b.
There is abundant seismicity on the eastern side of Coachella Valley at ≥6-km depth, located east of the surface trace of the San Andreas fault, which conforms with our inferred subsurface location of the San Andreas fault in Figure 8c. Seismicity patterns on the eastern side of Coachella Valley reveal a complex connection between the main San Andreas fault strand and east-west and north-south striking faults. At ≥6-km depth, the north striking Eureka Peak fault appears to extend in seismicity, at least, to near the Indio Hills fault. Similarly, the east-west striking Smoketree Wash fault extends in seismicity to near the subsurface location of the San Andreas fault (e.g., Figures 2 and 8c).

Estimated Basement Surface in Coachella Valley
The basement surface beneath Coachella Valley basin is not well known due to lack of publicly available, deep penetrating wells that reach the basement (e.g., Jennings & Hart, 1956 ;Proctor, 1968). Thus, we base our basement surface on observed seismic velocities of exposed basement rocks in the area (Fuis et al., 2017, see section 5.1).
Using the 4.50-km/s isovelocity contour as the basement surface, we infer that sediment thickness within the valley ranges from a maximum of~4 km at the southern end of the valley near the Salton Sea to less than 2 km at the northern end of the valley (Figure 11). At the locations of wells in the basin (Jennings & Hart, 1956;Proctor, 1968), which were all bottomed in either sand, gravel, or alluvium, our basement surface depth exceeds the maximum depth of the wells. In our inferred basement map, several features we discuss in the seismic velocities, such as Coachella Valley basin asymmetry, become evident. The basin geometry is also consistent with gravity derived basement depths from Langenheim et al. (2005) and Dorsey and Langenheim (2015) shown as yellow contours in Figure 2. However, the gravity-derived basement contours indicate sedimentary layer thickness up to 6 km. Two-dimensional gravity modeling near SSIP Line 4 from Dorsey and Langenheim (2015) estimate the maximum sedimentary thickness in the valley to be 4 km, but results from such models are known to be nonunique. Therefore, our results not only serve as a validation for previous gravity studies but are also consistent with previous SSIP studies.
The Peninsular Ranges block has a shallow dip beneath the valley. Fattaruso et al. (2014) and Dorsey and Langenheim (2015) indicate that the most probable explanation for the basin asymmetry is oblique convergence at the San Andreas fault, caused by differing orientations between the San Andreas fault strike and the Pacific-North American plate motion vector. This convergence acts as a vertical load, which in turn causes a rigid downward northeast tilting of the Peninsular Ranges block about a horizontal axis parallel to the San 10.1029/2018JB016260

Journal of Geophysical Research: Solid Earth
Jacinto fault (Dorsey & Langenheim, 2015). Boundary element modeling of a northeast dipping San Andreas fault supports the Peninsular Ranges block rotation and localized uplift in Mecca Hills (Fattaruso et al., 2014).
Our basement surface shows the prominent northwestward step-up (labeled in Figure 11a), which we interpret in the seismic velocities (Figure 7b). This basement step divides the Coachella Valley basin sediments at the junction where the San Andreas fault splits into the Banning and Mission Creek strands (Figure 11). Other basement step-ups are apparent to the north (Figure 11b) and suggest that Coachella Valley may consist of multiple subbasins (Persaud et al., 2014). Finally, our basement surface in the mountain ranges ( Figure 11) has a good correlation with the surface topography ( Figure 2) and aligns well with surface traces of mapped faults in the area, especially the San Andreas fault.

Conclusions
We inverted local earthquakes and SSIP explosion travel time data in Coachella Valley to obtain a detailed 3-D velocity model. Our velocity model agrees with the surface geology and previous gravity and seismic studies. Lateral velocity changes across the main faults exist to 9-km depth, indicating that major faults penetrate deep into the crust. Seismic velocities in our model define Coachella Valley basin asymmetry and crustal heterogeneity in a much greater detail than previous local and regional 3-D tomographic models of the northern Salton Trough of comparable scale.
The basement in Coachella Valley is complex and is overlain by sediments that reach an estimated thickness of~4 km near Mecca Hills at the northern end of the Salton Sea, gradually decreasing to a thickness of~2 km at the northern end of the valley. Basin geometry is asymmetric with more sediment accumulation on the eastern side of the valley against the San Andreas fault, consistent with previous gravity and 2-D seismic refraction studies. Combined interpretation of seismicity, lateral velocity changes, and observation of geologic offset reveal a zone in the Little San Bernardino Mountains that can be interpreted as at least one previously unmapped fault striking northeast from the Eureka Peak fault, which we refer to as the Lost Horse Valley fault zone.
Our model has a practical significance for improving earthquake hazard studies by providing a more accurate seismic velocity model for the northern Salton Trough. The accuracy of ground motion estimates strongly depends on the seismic velocity structure, especially the basin structure which is key in determining shaking intensity . The Coachella Valley basin asymmetry and the irregular basement structure defined in our velocity model will result in different ground shaking estimates than for a symmetric basin and regular basement structure in current regional community velocity models used in seismic hazard analysis for Southern California. Strong crustal heterogeneities present in our model can be used to improve modeling of the rheological properties of the crust and their interaction with main faults in the area.