Western US intraplate deformation controlled by the complex lithospheric structure

Domain: earth_science
Published: 2024-05-09
Editorial outcome: major_revision

Abstract
The western United States is one of Earth’s most tectonically active regions, characterized by extensive crustal deformation through intraplate earthquakes and geodetic motion. Such intracontinental deformation is usually ascribed to plate boundary forces, lithospheric body forces, and/or viscous drag from mantle flow. However, their relative importance in driving crustal deformation remains controversial due to inconsistent assumptions on crustal and mantle structures in prior estimations. Here, we utilize a fully dynamic three-dimensional modeling framework with data assimilation to simultaneously compute lithospheric and convective mantle dynamics within the western United States. This approach allows for quantitative estimations of crustal deformation while accounting for the realistic three-dimensional lithospheric structure. Our results show the critical role of the complex lithospheric structure in governing intraplate deformation. Particularly, the interaction between the asthenospheric flow and lithospheric thickness step along the eastern boundary of the Basin and Range represents a key driving mechanism for localized crustal deformation and seismicity.

Introduction
The driving mechanism behind intraplate deformation remains a fundamental scientific question. One example is the western United States (WUS), which displays a broad and complex crustal deformation pattern (Fig. 1). Since the mid-Cenozoic, much of the WUS has been experiencing long-lasting crustal extension, with well-documented geological and geophysical evidence1. However, the forces driving the extensional tectonics remain debated, with proposed mechanisms broadly falling into three categories: long-range plate boundary forces2–4, lithospheric body forces resulting from lateral gradients of gravitational potential energy (GPE)2–8, and basal tractions exerted by the underlying mantle flow5,7,9–13. Many researchers realized that a single mechanism is insufficient to explain the observed crustal deformation1–5,7,9–11,14, but the importance of different mechanisms is still controversial. One key outstanding question is how the three-dimensional (3D) lithospheric structure, a feature usually simplified in previous studies, is linked to the proposed driving forces.

Contemporary intraplate deformation in the WUS is evident through the observed crustal motion deviating from rigid plate rotation and widespread intraplate earthquakes. In the WUS, Global Positioning System (GPS) measurements revealed a clockwise rotation pattern15 in lateral crustal motion relative to the stable North American plate16 (Fig. 1a). Generally, the Basin and Range (B&R) is continuously extending in the northwest direction17,18, while the Great Valley (GV) and Colorado Plateau (CP) remain relatively undeformed17–20. Notably, the crustal motion within the southern B&R is lower than that further north21. From the northwestern B&R to the Pacific Northwest (PNW), the crustal motion transitions from northwest to northeast15,22. This change in the surface motion was primarily attributed to the changing relative driving force, transitioning from right-lateral shear along the Pacific-North American plate boundary to northeastern push from subduction of the Juan de Fuca (JdF) plate15. While this interpretation aligns with observed surface kinematics, the broader deformation over the entire WUS may also reflect other driving forces, such as those from lateral gradients of lithospheric GPE and active mantle flow, whose exact roles remain unclear.

The seismicity in the WUS correlates with the geodetically measured crustal motion (Fig. 1a) and can be categorized into four main belts/regions: the San Andreas Fault (SAF), Walker Lane (WL), Intermountain Seismic Belt (ISB), and PNW. Recent geophysical studies revealed that the widespread seismicity in the WUS is strongly associated with the spatially varying crustal stress state23–26 (Fig. 1b). The relatively undeformed GV is bounded by the SAF and WL with prominent seismicity, characterized mainly by strike-slip faulting (Fig. 1b). The remarkable ISB, a north-south-trending seismic belt stretching from the southern CP to the northern Rocky Mountain (RM)27,28, coincides with the eastern boundary of B&R as well as the transition from thin to thick North American lithosphere (Fig. 1a). Most earthquakes in the ISB show normal faulting with E-W extension along the eastern boundary of B&R and the tectonic parabola surrounding the Yellowstone hotspot track, but the crustal stress transitions to a more compressional state in the northern RM with earthquakes exhibiting strike-slip to thrust faulting (Fig. 1b). Seismicity in the PNW is spatially distributed, mainly in the back-arc region, which corresponds to the diffuse crustal motion in this area.

The debated mechanisms for WUS crustal deformation could be boiled down to the often simplified and yet uncertain lithospheric and convective mantle structures adopted in different studies and the associated fine-scale dynamics2–14,29–37. Particularly, the lithosphere–asthenosphere interaction below this region remains contentious. Previous studies generally assumed a flat lithosphere–asthenosphere boundary (LAB) for the continental plate at 100 km or deeper2–7,9–11,13. Under this assumption, the convective mantle and lithospheric dynamics were often computed separately and connected through one-way coupling that only allowed the convecting mantle to exert horizontal shear (i.e., viscous drag) along the flat LAB4,6,7,11,12. This horizontal shear was usually estimated from mantle convection models that only considered long-wavelength deep density anomalies4,6,7,11,12. Consequently, this approach may not adequately resolve the fine-scale interaction between the convective mantle and the lithosphere due to the lack of spatially varying lithospheric thickness and associated mantle deformation35–37. Indeed, the assumed flat LAB contradicts recent global38,39 and regional40–42 studies showing significant lateral variations in lithospheric thickness.

In addition, the importance of 3D lithospheric effective viscosity structure, which impacts stress transmission and crustal deformation, has been previously undervalued. Earlier studies either assumed mechanically strong continental plates with uniform strength9,11–13 or focused solely on lateral variations in bulk lithospheric effective viscosity2–5,7,10–12. In the latter scenario, lithospheric dynamics was frequently modeled under the thin-shell approximation, treating the lithosphere as a thin two-dimensional (2D) layer with constant thickness and vertically uniform properties. This method usually estimated lateral viscosity variations from geodetically or geologically determined deformation fields2,3,5,7,10. However, recent studies highlighted the crucial need for independently determined 3D lithospheric effective viscosity structures to accurately model crustal deformation31. Hence, a more sophisticated modeling approach beyond the thin-shell approximation is essential to fully capture the detailed 3D lithospheric dynamics.

Results
Construction of the fully dynamic 3D model
We model the lithosphere and convective mantle as an incompressible viscous fluid in spherical geometry (“Methods” section). The numerical model covers a region spanning 60°×100°×2890 km in latitude×longitude×radius, with the finest resolution (~7 km) inside the lithosphere and asthenosphere within the WUS. The model’s high resolution allows us to faithfully simulate the fine-scale 3D dynamics in the WUS, especially the previously unexplored lithosphere–asthenosphere interaction arising from the spatially varying lithospheric thickness and shallow mantle flow. Unlike most previous studies that employed one-way coupling, our model incorporates natural coupling between the convecting mantle and lithosphere by simultaneously including the realistic lithospheric and convecting mantle structures (Fig. 2). Furthermore, the data-oriented nature of our geodynamic modeling approach allows us to test models with different lithospheric and convecting mantle structures to isolate their respective effects on crustal deformation. In this study, we construct ten different models (Table 1) to fully explore the driving mechanism for WUS intracontinental deformation.

In our model (“Methods” section), the present-day mantle structure is adopted from a time-dependent mantle flow model using the hybrid geodynamic modeling approach44, which combines the forward45 and adjoint46 data-assimilation algorithms. The present-day 3D JdF slab structure was reproduced by the high-resolution forward approach45, and other mantle structures, including the subducted Farallon slab, lithospheric drips, and hot-mantle anomalies, were derived from recent high-resolution seismic tomography42,47 via the adjoint method46. This time dependence matches several independent geophysical and geological observations within the WUS. For example, the evolving mantle thermal state reproduced the late-Cenozoic volcanic history36 and topographic changes48. The associated mantle flow and density anomalies also matched the SKS splitting measurements37 and free-air gravity anomaly48.

The density and viscosity structures of the continental lithosphere, characterized by strong lateral variations in thickness, play a critical role in truthfully modeling lithospheric deformation2–4,6,7,14,31. We employ a hybrid approach that merges LITHO 1.038 and a recent regional body-wave tomography model42 to construct our lithospheric thickness model (“Methods” section, LTM1), which effectively captures the first-order lithospheric thickness change along the eastern boundary of B&R (Fig. 1). Using the seismically detected Moho depth49 and LAB depth38,42, we deduce the lithospheric density structure by matching surface topography (“Methods” section). Additionally, we construct a 3D lithospheric effective viscosity structure based on recent seismic attenuation studies50,51 (“Methods” section, Fig. 2g). This approach allows us to independently calculate the fine-scale lithospheric deformation by avoiding direct reliance on GPS measurements.

Lithosphere–asthenosphere interaction beneath the WUS
Our fully dynamic 3D model well captures the fine-scale interaction between the lithosphere and convecting mantle beneath the WUS. The dominantly eastward upper-mantle flow in the region is mainly driven by the sinking Farallon slab currently located beneath the central and eastern US36,37 (Fig. 2a–f). Locally, the flow is diverted by the subducting JdF slab and the thick cratonic lithosphere on the east (Fig. 3a). Below the region with a thin lithosphere, the predicted mantle flow pattern independently matches the recent observation of asthenospheric azimuthal anisotropy52 (Fig. 3a). The modeled highly segmented JdF slab, consistent with seismic tomography53, allows hot asthenospheric material to flow through the central slab tear and around its edges to beneath the northern B&R and southern SRP36,37, respectively. Further inland, the eastward flow is bisected by the thick CP lithosphere, where the southern branch of hot mantle flows around the CP into the Rio Grande Rift (RGR), and other hot material flows northeastward until blocked by the thick cratonic lithosphere east of the B&R. This peculiar flow pattern is vividly tracking the azimuthal anisotropy below the thin WUS lithosphere, with an average angular misfit of 37.5° between the directions of flow and anisotropy (“Methods” section). In the cratonic region, our model predicts little internal deformation within the lithospheric mantle. This implies that there is no active shearing within the lithosphere to generate prominent azimuthal anisotropy. Therefore, we interpret the sudden change of anisotropy orientation across the eastern boundary of B&R as reflecting a contrast between active asthenospheric deformation on the west and frozen-in fabric in the cratonic lithospheric mantle to the east.

To illustrate the importance of lateral variations in lithospheric thickness on WUS mantle dynamics, we analyze the asthenospheric flow using different lithospheric thickness models (Fig. 3b). The predicted asthenospheric flow exhibits near-identical patterns with a laterally varying (LTM1; red arrows in Fig. 3b) and a constant (LTM2; light-blue arrows in Fig. 3b) lithospheric thickness beneath the region with a thin lithosphere. However, the fast horizontal flow with LTM2 continues to go eastward beneath the intermountain west, while the flow in the model with LTM1 largely stops here (Fig. 3b). Notably, the fast eastward mantle flow predicted with LTM2 is nearly orthogonal to the observed azimuthal anisotropy beneath the intermountain west, resulting in a large average angular misfit of 49.7°. Although seismic anisotropy does not uniquely reveal the pattern of mantle flow, it could help eliminate certain scenarios. For example, the flow predicted with LTM2 (Fig. 3b) suggests strong E-W azimuthal anisotropy east of the B&R, opposite to observation (Fig. 3a). In contrast, the flow with LTM1 implies little deformation in this region, not violating the anisotropy constraint that could reflect fossil fabric instead.

Complex lithospheric structure controls WUS crustal deformation
The E-W contrast in WUS lithospheric thickness (Fig. 1) implies a complex 3D effective viscosity and density structure. In the thin lithosphere region on the west, both the crust and lithospheric mantle are rheologically weak (Fig. 2d–g) and actively deforming (Fig. 1a). In contrast, the cratonic lithosphere east of B&R is much thicker and less tectonically active; this thick root blocks and diverts shallow asthenospheric flow (Fig. 3a). Besides, recent studies demonstrated that the thick cratonic lithospheric mantle is denser than the ambient asthenosphere and modulates the crustal stress and surface topography14,54–56. Despite these findings, the role of 3D lithospheric structure in modulating crustal deformation remains largely unexplored.

To quantitatively assess the impact of 3D lithospheric structure on crustal deformation, we conduct three simulations (M1-M3) using different lithospheric structures with increasing complexity. The crustal deformation rate (i.e., the second invariant of strain rate tensor) and seismicity are expected to be correlated57–59. Therefore, we compare the predicted crustal deformation rate to observed seismicity, while employing a seismicity prediction power analysis (“Methods” section) to quantify the match between model prediction and observation, with a higher value representing a better match. To compare with the GPS measurements, we calculate the vectorial mismatch between the predicted surface motion and a smoothed version of GPS-determined surface motion (“Methods” section).

The model with the simplest lithospheric structure (M1), consisting of a laterally uniform crust and lithospheric mantle, predicts a diffusive deformation pattern without clear localization of strain rate (Fig. 4a). The seismicity prediction power of M1 is 34.33%. Most modeled deformation is concentrated along the plate boundary (i.e., the SAF). The predicted crustal motion, whose magnitude decays inland (Fig. 4d), is dominated by the plate boundary effect (Supplementary Fig. 4a–c). The predicted surface motion deviates from observation by an average residual of 4.01 mm/yr.

In the second model (M2, Fig. 4b, e), we further consider the lateral gradients of crustal GPE, another primary driving force for intraplate deformation3,8,60. We approximate this GPE using laterally varying Moho depth60 and a constant crustal density at 2850 kg/km3, except for the SRP, whose density is assumed to be 2950 kg/m3 due to its enrichment in basaltic composition61,62. Other model properties in M2 are the same as in M1. The resulting crustal deformation rate is slightly elevated in regions with thick crust, such as the CP and Wyoming Plateau (WP), due to locally elevated GPE (Fig. 4b). Additionally, the high GPE in the intermountain west drives westward crustal motion in the WUS (cf. Fig. 4d, e). M2 has a seismicity prediction power of 35.50% and an average residual in crustal motion of 4.44 mm/yr, both similar to M1. Clearly, neither M1 nor M2 reproduces the observed seismicity distribution with localized strain accumulation and crustal motion. These tests call for a heterogeneous lithosphere, likely with laterally varying thickness and 3D variations in density and effective viscosity structures14,31,54–56, to fully reproduce the observed crustal deformation field.

By incorporating the geophysically inferred lithospheric structure (Figs. 1 and 2), the resulting crustal deformation pattern correlates much better with the observed seismicity distribution (M3, Fig. 4c) and crustal motion (Fig. 4f), with the seismicity prediction power dramatically increasing to 96.37% and the average crustal motion residual dropping to 1.78 mm/yr. In particular, the WL and ISB display very high deformation rates, consistent with their intense seismicity. In contrast, the GV and CP have low deformation rates, also agreeing with seismicity and geodetic measurements17–20. The high deformation rate in the PNW extends far inland to the western SRP, supported by the distributed seismicity in this region. In addition, the observed rotational crustal motion pattern in the WUS is well reproduced in this model (Fig. 4f). Therefore, the heterogeneous lithosphere is important for forming the observed crustal deformation pattern.

Estimations of different driving forces and associated deformation
The estimations of proposed driving forces, namely plate boundary forces, lithospheric body forces, and basal tractions, depend on the adopted lithospheric structures. In our model, the plate boundary effect is defined as forces exerted on the continental plate by the moving oceanic plates with no contribution from the underlying mantle flow. Traditionally, basal tractions were defined as the horizontal shear exerted by mantle flow along a flat LAB4,6,7,11,12. However, in reality, the LAB exhibits undulations (Fig. 1) and dynamically interacts with the convecting mantle (Fig. 3a), resulting in complex surface tractions comprising both shear and normal forces4. Thus, in our model, we define basal tractions as the total surface forces along the undulating LAB, originating from the underlying mantle flow.

The conventional estimation of lithospheric body forces relied on computing lateral gradients of lithospheric GPE, assuming lithospheric isostasy with surface topography fully supported by lithospheric density anomalies above a presumed LAB depth of 100 km2,3,7. However, this assumed global compensation depth might not be appropriate, as the convective mantle also contributes to surface topography via radial stress along the LAB (i.e., dynamic topography). Furthermore, recent seismic studies revealed strong lateral gradients in LAB depth38–42, particularly in regions like the CP and WP (Fig. 1), where the lithosphere exceeds 100 km in thickness. Consequently, the traditional calculation of lithospheric GPE is physically inaccurate. In our fully dynamic 3D model, lithospheric body forces refer to those originating from lithospheric density anomalies down to the undulating LAB. Unlike previous studies that take surface topography as an input for lithospheric GPE calculation2,3,5,7–12, we use it as an independent observation to validate our lithospheric density structure (“Methods” section).

We construct seven additional geodynamic models (M4–M10, illustrated later) with different lithospheric and convective mantle structures (Table 1) to fully explore the impact of the heterogeneous lithosphere on crustal deformation and to quantitatively analyze the different driving forces. In all models, plate motion63 is prescribed in the ocean basins, stable central-eastern US and Canada, as velocity boundary conditions (“Methods” section). Hence, the predicted deformation rate and surface motion include the plate boundary effect in all models.

Mantle convection influences lithospheric deformation through two mechanisms: the plate boundary effect, reflecting the integrated effect of global mantle flow, and basal tractions, arising from mantle flow right beneath the plate. To evaluate how a heterogeneous lithosphere impacts the contribution of mantle convection to crustal deformation, we perform six simulations and analyze the resulting crustal deformation patterns (Fig. 5 and Supplementary Fig. 5). In the models solely considering the plate boundary effect (M4, M7, and M8), we exclude all density anomalies in the convective mantle and continent lithosphere (Fig. 5a, d, Supplementary Fig. 5a, b, e, f). To further consider basal tractions, we utilize three models (M5, M9, and M10), each incorporating different lithospheric structures and all density anomalies in the convecting mantle (Fig. 5b, e, Supplementary Fig. 5c, d, g, h). The crustal stress shown in Fig. 5b reflects the difference between M5 and M4, which have the same lithospheric structure but different convective mantle structures, thereby highlighting the contribution of basal tractions to crustal stress.

We first estimate the crustal stress and deformation resulting from the plate boundary effect by constructing M4, which assimilates a geophysically inferred heterogeneous lithosphere consisting of LTM1 and 3D variations in effective viscosity. Compared to the models solely considering the plate boundary effect but having a constant lithospheric thickness (M7 and M8, Supplementary Fig. 5a, b), M4 has a larger seismicity prediction power (89.88%). Notably, along the lithospheric thickness step (the orange dashed line in Fig. 5a), the predicted crustal deformation is localized on the thin and weak side, with little crustal deformation in the interior of CP and WP (cf. Fig. 5a, Supplementary Fig. 5a, b). This finding highlights the importance of the lithospheric thickness step in generating localized crustal deformation. East of the lithospheric thickness step, the thick and strong lithospheric keel penetrates deeper into the asthenosphere, thereby enhancing the coupling between the thick lithosphere and the underlying static mantle. This locally enhanced coupling slows the northwestward motion of the thick lithospheric root (cf. Fig. 5d, Supplementary Fig. 5e, f, magenta arrows), effectively counteracting the shearing from the Pacific plate. As a result, the average residual reduces to 1.93 mm/yr.

In model M5, which further incorporates basal tractions, the predicted crustal deformation exhibits clearly localized high deformation rates in some intraplate regions, such as in the WL and ISB (Fig. 5b), with a seismicity prediction power of 96.62%. The inclusion of shallow hot anomalies (between the LAB and 100 km depth) that are not considered in cases incorporating basal tractions but having a flat LAB at 100 km (M9 and M10) leads to a reduction in asthenospheric effective viscosity (Fig. 2g). Consequently, the coupling between the lithosphere and convective mantle is locally weakened, resulting in a reduction in deformation rate in the interior of northern and southern B&R (cf. Fig. 5b, Supplementary Fig. 5c, d). In contrast, along the lithospheric thickness step (the orange dashed line in Fig. 5b), the crustal deformation rate is largely unchanged (cf. Fig. 5b and Supplementary Fig. 5d), indicating the interaction between the mantle flow and thick lithospheric keel generates localized crustal deformation along the step that counteract the decoupling effect. Compared to the case without the density-driven mantle flow (M4, Fig. 5a), it shows that the interaction between the active mantle flow and thick lithospheric root locally increases the crustal deformation rate, like in the northern RM and along the eastern boundary of northern B&R (cf. Fig. 5a, b). The eastward mantle flow also pushes the lithospheric root eastwards, reducing westward motion in the southern CP and central and southern B&R (cf. Fig. 5e, Supplementary Fig. 5g, h).

To estimate crustal stress and deformation due to the lithospheric body forces over the WUS, we construct a model (M6) that incorporates a heterogeneous lithosphere with 3D viscosity and density variations. In this model, all density anomalies in the convective mantle are removed. The crustal stress shown in Fig. 5c is the difference between M6 and M4, representing the crustal stress solely due to the lithospheric body forces. The lithospheric body forces, arising from lateral gradients of lithospheric GPE, elevate the crustal deformation rate over the whole WUS (cf. Fig. 5a, c). While our estimated crustal stress due to the lithospheric body forces is generally consistent with previous studies2–4,7–9, it is noteworthy that the underlying physics differs. We emphasize that a continental lithosphere with a laterally varying thickness and a denser-than-asthenosphere lithospheric mantle is critical in calculating and understanding the lithospheric GPE14,56. Additionally, the high GPE east of the B&R drives the WUS to move to the northwest, particularly in the B&R, SRP, and PNW (cf. Fig. 5d, f).

Mechanism for localized intraplate deformation
While the 3D crustal effective viscosity structure influences the crustal deformation rate (“Methods” section), the primary driving force for the localized crustal deformation along the lithospheric thickness step is the locally enhanced interaction between the lithosphere and asthenosphere (Figs. 5 and 6). This enhanced lithosphere–asthenosphere interaction ultimately affects the crustal effective viscosity structure, particularly the intraplate crustal weak zone (Fig. 2g), which cannot be fully explained by other mechanisms (“Methods” section). Our model demonstrates that beneath regions with a thin lithosphere, the mantle flow does not generate significant shear on the horizontal plane (Fig. 6a–c) due to the low asthenospheric effective viscosity (Fig. 2g), which is consistent with previous studies6,7. However, our model reveals a strong eastward normal traction (i.e., pressure drag) on N-S trending vertical planes beneath the region with a thin lithosphere (Fig. 6e, f), which is not previously captured. This increased eastward normal traction is due to the thick lithospheric keel blocking the landward mantle flow (M3, Figs. 3 and 6), as it leads to large positive dynamic pressure in the convecting mantle beneath the thin lithosphere region. This pressure drag (~20 MPa), which is part of basal tractions, pushes the lithospheric keel eastwards and increases the horizontal shear in regions with a thick lithosphere, such as the CP, WP, and northern RM (cf. Fig. 6a–c). Furthermore, this pressure drag drives localized lithospheric deformation along the lithospheric thickness step at the eastern boundary of B&R (Fig. 5a, b), representing a key mechanism for the ISB (Figs. 1a and 4).

Recent seismological studies41 revealed fine-scale lithospheric thinning beneath the eastern boundary of northern B&R, a fact that is not included in our lithospheric thickness model. This observed lithospheric thinning may lead to further localization of lithospheric deformation by reducing the bulk lithospheric strength, similar to the effect of a localized crustal weak zone. The shallow LAB observed beneath the eastern boundary of northern B&R appears anti-correlated with the crustal thickness, suggesting its formation as a result of recent tectonic activities41. Therefore, we propose that the localized interaction between hot asthenospheric flow and lithospheric thickness step, starting around 12 Ma36, plays an important role in generating locally thinned lithosphere.

Revised role of different forces in driving WUS crustal deformation
The importance of variable lithospheric properties (Figs. 3–6) requires a reconsideration of previously proposed driving forces (plate boundary forces, lithospheric body forces, and basal tractions) for WUS crustal deformation. In many previous studies, the estimations of driving forces were influenced by the simplified lithospheric structures adopted2,3,5–7,10–12. We have shown that the lithospheric thickness step creates a strong pressure drag (Fig. 6) that causes localized intraplate deformation (Figs. 4 and 5), a factor that has been previously neglected. Additionally, the presence of a dense and thick cratonic root alters the estimation of lithospheric GPE, which calls for a reevaluation of its contribution to crustal deformation. Our best-fitting model (M3) not only successfully reproduces the observed crustal deformation (Fig. 4c, f), but also predicts a crustal stress field consistent with high-quality earthquake focal mechanism solutions23,24 and recent geophysical estimations of the horizontal maximum compression directions (SHmax) and Aφ parameter25,26 (Fig. 7). Based on the best-fitting model, we further evaluate the relative importance of plate boundary forces (Fig. 5a), basal tractions (Fig. 5b), and lithospheric body forces (Fig. 5c) by calculating the ratios of their magnitudes (i.e., the second invariant of stress tensor; Fig. 8).

The observed focal mechanism solutions23,24, geophysically estimated SHmax directions25,26, and Aφ parameter25,26 show strong spatial variations over the WUS that cannot be fully explained by one single driving force9,14,25,26. In the southwest coastal region, the plate boundary effect dominates the crustal motion and stress field (Figs. 4f, 5a, d, 7a and 8e, f). Along the SAF, the Pacific plate exerts right-lateral shearing forces, driving strike-slip faulting (Figs. 5a and 7a, c) and northwestward crustal motion (Fig. 5d). The PNW and northern SRP are complex regions with observed crustal stress generally in N-S compression23–26 (Fig. 1b), which was largely attributed to the plate boundary forces from the shearing along the SAF and subduction of the JdF plate15,25,64. In our best-fitting model, all three driving mechanisms contribute large stresses in the PNW (Fig. 5a–c). The crustal stress generated by GPE gradients and basal tractions is generally in the same direction with opposite signs in this region. Combined with the plate boundary forces exerted by the JdF plate and Pacific plate, the PNW and northern SRP experiences strong N-S to NW-SE compression (~20–40 MPa), with crustal stress state in the thrust faulting regime (Fig. 7c). Overall, the basal tractions play the most important role in deforming the PNW and northern SRP (Fig. 8d–f). Further inland, in the central and southern B&R, the westward extension (Fig. 7a, c) is predominantly driven by the plate boundary force originating from the Pacific-North American boundary with strong E-W extension and weak N-S compression (Fig. 5a). Furthermore, the N-S extensional stress due to basal tractions (Fig. 5b) and E-W extensional stress from the gradients of lithospheric GPE (Fig. 5c) collectively modify the crustal stress field into purely E-W extension, aligning with the focal mechanism solutions and geophysical estimations of SHmax directions and Aφ parameter (Fig. 7). In the northern B&R, extensional crustal stresses due to plate boundary forces and GPE gradients generally have comparable magnitudes (Fig. 8f). In the northwestern B&R, the crustal stress field shows a combined effect of all three driving forces with comparable magnitudes (Fig. 8d–f). Further north, the E-W compression in the northern RM is predominantly from basal tractions (Fig. 8d, e), which are locally enhanced by the strong and neutrally buoyant slab curtain65. Noteworthy, along the eastern boundary of B&R and the tectonic parabola surrounding the Yellowstone hotspot track, where the ISB is located, the GPE-induced extensional stress dominates the crustal stress field and drives normal to strike-slip faulting (Figs. 7a and c, 8d, f), in line with the focal mechanism solutions23,24. Overall, in the WUS, the plate boundary effect controls the coastal region and the central-southern B&R, while lithospheric body forces dominate the northern B&R, and basal tractions resulting from mantle flow play an important role in deforming the northern WUS (Figs. 8 and 9).

Our systematic analysis demonstrates the crucial role of a 3D heterogeneous lithosphere fully coupled with the underlying convecting mantle in accurately reproducing observed crustal deformation. In particular, the lateral variation in lithospheric thickness significantly influences the contribution of mantle convection to crustal deformation, as well as the estimation and understanding of lithospheric body forces (Figs. 5, 9, and Supplementary Fig. 5). The GPS-measured lateral surface deformation in the WUS reflects a combined effect of plate boundary forces, lithospheric body forces, and basal tractions. Among all driving forces, only basal tractions drive eastward differential motion over the WUS, a crucial factor in reproducing the observed surface deformation (Figs. 4f and 5d–f). This finding highlights the importance of shallow mantle flow interacting with a lithosphere that has laterally varying thickness (Figs. 3 and 4), as is often neglected in previous studies. Furthermore, the heterogeneous lithospheric structure plays a critical role in generating the observed crustal deformation within the WUS. Particularly, the lithospheric thickness step along the eastern boundary of B&R is essential in generating localized intraplate deformation, which leads to the formation of the enigmatic ISB (Figs. 1 and 6).

Methods
Data-oriented geodynamic modeling approach
To construct a self-consistent, fully dynamic model with coupled lithosphere and convective mantle, we generally follow the modeling approach described in Cao and Liu14. The lithosphere and mantle are modeled as an incompressible viscous fluid in spherical geometry under the Boussinesq approximation. We use a user-modified version of CitcomS with tracers14,66,67 to solve the momentum equation. The high-resolution 3D models cover a region of 60°×100°×2890 km in latitude×longitude×radius, with the finest resolution (~7 km) inside the lithosphere and asthenosphere within the WUS. This model domain is sufficiently larger than the conterminous US. It contains the surrounding ocean basins and parts of Canada and Mexico, allowing natural modeling of the detailed 3D lithospheric and mantle dynamics, including the fine-scale lithosphere–asthenosphere interaction that was not captured in previous studies.

A hybrid velocity boundary condition is used on the top of the model to incorporate far-field forces and truthfully simulate intracontinental deformation. Plate motion63 is prescribed in the ocean basins surrounding the conterminous US, the stable central and eastern US (east of 110°E), and part of Canada (north of 55°N), while the actively deforming WUS and Mexico are set to have free-slip boundary condition. All other boundaries have a free-slip velocity boundary condition. Each model is an individual simulation with different lithospheric and mantle structures, as listed in Table 1.

Lithospheric and convective mantle structures
The convective mantle structure is adopted from a time-dependent mantle flow model using the hybrid geodynamic modeling approach44. This time-dependent convecting mantle structure matches multiple geophysical and geological observations36,37,48. Most features in the convecting mantle, including the subducted Farallon slab, lithospheric drips, and hot asthenospheric material, were derived from high-resolution seismic tomography models42,47 using the adjoint assimilation scheme46. Features near the subduction zone, including the subducting Juan de Fuca slab and mantle wedge, were generated by a forward sequential assimilation scheme45. The convecting mantle has a temperature- and depth-dependent effective viscosity structure, the same as Zhou et al.36,37.

The lithosphere structure in our model is inferred from recent seismic studies38,42,49–51,68–70. The initial lithospheric thickness model has a Moho depth model from CRUST 1.049 and a hybrid LAB depth model14 derived from LITHO 1.038 and a body-wave tomography model42. Since LITHO 1.038 does not capture the abrupt lithospheric thickness change along the eastern boundary of the B&R, which is better imaged in the recent regional body-wave seismic tomography42, we redefine this part of the LAB depth based on the body-wave seismic tomography model. In practice, we assume the continuous fast anomalies down to 250 km from the seismic tomography model represent the continental lithosphere. We also adopt an MLD in the thick cratonic lithospheric mantle at 120 km, based on recent seismic studies68–70. Seismic studies revealed that the MLD marks the boundary between two layers with different seismic properties68–70, also indicating changes in their physical properties. Besides, recent geophysical and geodynamic studies showed that the continental lithospheric mantle has a layered density structure14,54–56,71–73. Therefore, we assume the MLD also marks an interface for lithospheric mantle densities.

The initial lithospheric density structure has three layers. In the crust, we assume a constant density of 2850 kg/m3, except for the SRP. In the SRP, its lower-than-ambient elevation was attributed to its basalt-rich crust61,62, so we adopt a crustal density of the SRP at 2950 kg/m3, which is the same as that in Zhou and Liu48. The density of the lithospheric mantle is elusive, and the proposed density ranges from neutrally buoyant74–76 to denser than ambient asthenospheric material14,54–56,71–73,77. In our lithosphere model, the density structure of the lithospheric mantle has two layers separated by the seismically observed MLD68–70. The densities of these two layers were adjusted to reproduce the continental-scale surface topography contrast between the WUS and the central and eastern US, as well as between the Pacific Ocean and the WUS14. The resultant best-fitting lithospheric mantle density structure has an upper layer with a density of 3355 kg/m3 over the whole conterminous US and a lower layer with a density of 3375 kg/m3 in the central US and 3385 kg/m3 in the eastern US, respectively. This lithospheric mantle density structure is consistent with recent geophysical estimations56,72,73, which showed the cratonic lithospheric mantle is overall ~1% denser than ambient asthenosphere (3340 kg/m3).

To capture the fine-scale lithospheric GPE variations, we fix the density for each layer and region, then adjust the LAB and Moho depth to reproduce the surface topography within the WUS. This process involves two steps: (1) calculate the initial topography difference between the observation and prediction, then adjust the LAB depth to accommodate the initial topography difference with constraints detailed in Cao and Liu14; (2) calculate the residual topography difference between the observation and new prediction from step 1, then adjust the Moho depth to accommodate the residual topography difference. In each step, we first conduct a simulation with the unmodified Moho and LAB depth to predict the surface topography. In the simulation, all density anomalies in the convecting mantle are included, so dynamic topography is naturally included in the prediction. This modeling approach does not rely on the assumption of isostasy. The resulting Moho depth is generally consistent with the recent seismically determined Moho depth in the WUS49,78 (Supplementary Fig. 6).

The 3D effective lithospheric viscosity structure is critical in modeling lithospheric dynamics3,4,6,7,14,31. In our lithosphere model, we use a modified version of the 3D effective lithospheric viscosity structure from Cao and Liu14. The continental-scale lithospheric viscosity structure was determined by a parameter search14. In the parameter search, the effective viscosity of the WUS, central US, and eastern US were independently varied. The best-fitting lithospheric effective viscosity structure minimized the angular misfit of maximum horizontal compression (SHmax) directions to the observation25,26 over the conterminous US.

We further infer the fine-scale lithospheric effective viscosity structure in the WUS based on recent seismic attenuation studies50,51. Seismic attenuation represents an elastic properties of rocks79,80 and could be linked to the steady state viscosity by linear viscoelastic theory. The same set of state variables controls both effective viscosity and seismic attenuation81,82. Generally, the seismic attenuation map determines the spatial pattern of effective viscosity but leaves the absolute value undetermined. Therefore, it is possible to infer the effective viscosity structure from seismic attenuation follows the same procedure as in Liu and Hasterok31. Since the seismic attenuation in the lithospheric mantle beneath the conterminous US shows regional-scale variations51, largely consistent with the effective viscosity variations in Cao and Liu14, we keep the effective viscosity in the lithospheric mantle unchanged. The fine-scale 3D crustal viscosity structure is converted from a recent crustal seismic attenuation study50. The crustal seismic attenuation shows the integrated effect through the whole crust. Practically, to construct the 3D effective viscosity structure, we first assume a linear crustal geotherm with lateral varying crustal thickness within the WUS. Then, we determine the temperature-dependent effective viscosity as a reference viscosity structure. Finally, we use the crustal attenuation map to introduce lateral variations at all depths, assuming seismic attenuation is evenly distributed in the radial direction. The whole process could be described by the following equations:

η = β·η0·exp(Ec·(T0+Tc)/(T+Tc))

β = a·(Q0/Q)

where η is the effective viscosity, β is a coefficient controlling the magnitude of lateral variations, η0 is the reference viscosity, Ec is the non-dimensional activation energy, T is the non-dimensional temperature, T0 is the surface temperature, Tc is the characteristic crustal temperature, Q is the crustal attenuation, Q0 is the reference crustal attenuation, a and c are numerical coefficients. Tc, η0, and Q0 are fixed in our parameter search. An extensive parameter search is performed to find the best-fitting conversion coefficients that minimize the averaged angular misfit to the observed SHmax directions25,26 (Supplementary Fig. 7). Through the test, we find the best-fitting coefficients are a=10, Ec=3, and c=3.

We want to emphasize that this approach only approximates the real 3D effective viscosity structure. A better way is to use the 3D magnetotelluric (MT) data inferred effective viscosity structure, which is critical in understanding lithospheric dynamics31,83,84. However, due to currently limited MT data coverage, it does not allow a survey covering the whole conterminous US. In the WUS, our effective viscosity model captures the main features revealed by a most recent MT survey85.

Calculating crustal stress and deformation
Our numerical models output the full stress tensor and the second invariant of the strain rate tensor on every computing node. We radially integrate the full stress tensor and the second invariant of the strain rate tensor down to 30 km (10 km for Supplementary Fig. 3). Subsequently, we calculate the depth-averaged horizontal principal stress and the crustal deformation rate (i.e., the second invariant of the strain rate tensor).

The predicted surface motion is output in a lower mantle fixed frame. We take the direct output and remove the North American plate motion63 relative to the fixed lower mantle to get the relative crustal motion.

Evaluating predicted mantle flow against seismic observation
Our convecting mantle structure is adopted from a time-dependent mantle flow model44, which matches multiple geophysical and geological observations36,37,46. After incorporating a geophysically constrained lithosphere with lateral varying thickness, we further evaluate the predicted mantle flow by comparing the active mantle flow to the recent observation of depth-dependent seismic azimuthal anisotropy52, which is not included in our previous studies. We define the active mantle flow as the asthenospheric material having a velocity larger than 0.01 cm/yr in the North American plate fixed reference frame16. Subsequently, we calculate the angular misfit between the predicted active mantle flow field and observed seismic azimuthal anisotropy52.

Seismic prediction power analysis
We perform a seismic prediction power analysis to compare the predicted crustal deformation rate field and observed seismicity. First, we calculate the normalized deformation rate field, which keeps the spatial pattern, by normalizing the prediction with the maximum predicted value within the WUS. Then, we sample the normalized crustal deformation rate at 26,577 earthquake epicenters during 1980–2020 (Fig. 1a). The quantitative relation between crustal deformation rate and seismicity remains unknown. Nevertheless, recent geodetic studies showed that more than 90% of seismicity happened in areas with geodetically measured crustal strain rate between 2×10−9 to 2×10−7 yr−1 (ref. 59). Therefore, we assume the model successfully predicts an earthquake if the normalized crustal deformation rate is larger than 0.01, within two orders of magnitude from the maximum, at the epicenter. Finally, we calculate the success percentage to quantify the match between the predicted crustal deformation rate field and observed seismicity.

Evaluating predicted crustal stress and surface motion
To quantify the match between the predicted crustal stress field and the observed crustal stress patterns, we calculate the SHmax directions and compare them to the observation25,26 following the same method as Cao and Liu14. We also calculate the Aφ parameter flowing ref. 86 to determine the faulting regime predicted by the model. Since the Aφ parameter is less well-constrained than SHmax directions26, we only qualitatively compare the prediction and observation (Fig. 7c).

The predicted crustal motion is compared to a smoothed version of GPS measurements (spline interpolated from the raw data) on a 0.5°×0.5° grid over the WUS (black arrows in Figs. 4d–f and 5d–f). The residual velocity, the vectorial difference between the prediction and observation, is calculated on each node, and its average magnitude and average angular misfit are used to quantify the match between the prediction and observation.

Effects of 3D crustal effective viscosity structure
We perform four additional models (M7–M10) with different lithospheric and convective mantle structures to quantify how the 3D crustal effective viscosity structure affects the contribution of mantle convection to crustal deformation. With a uniform lithosphere (M7), consisting of LTM2 and constant lithospheric effective viscosity, the predicted crustal deformation shows diffusive patterns and decay monotonically with the distance away from the plate boundaries (Supplementary Fig. 5a), leading to a seismicity prediction power of 73.10%. Within the regions far from the plate boundaries, such as the northern RM and WP, negligible deformation occurs. Besides, the predicted surface motion shows significant northwestward motion in the southwestern US and northward to northeastward motion in the PNW (Supplementary Fig. 5e), with a large average angular misfit (46.7°) and a significant average residual (4.50 mm/yr). The predicted clockwise rotation generally reflects the strong coupling with the fast-moving Pacific plate and the subducting JdF plate. Most of this clockwise motion is solid-body rotation, which does not involve internal deformation of the plate (cf. Supplementary Fig. 5a, e).

With the inclusion of basal tractions (M9), the predicted crustal deformation still maintains a smooth spatial pattern (Supplementary Fig. 5c). The seismicity prediction power increases to 78.11%. The subducting JdF slab locally enhances basal tractions along the western boundary of SRP, resulting in localized increments in crustal deformation (cf. Supplementary Fig. 5a, c). Overall, the active mantle flow drags the WUS eastwards and induces eastward surface motion, particularly in the PNW, SRP, and RM (cf. Supplementary Fig. 5e, g, magenta arrows). The active mantle flow reduces the angular misfit between predicted and observed surface motion to 34.8°, while the average residual remains similar (4.13 mm/yr in M9).

Subsequently, we introduce a 3D crustal effective viscosity structure into the model and assess its impact on the plate boundary effect (M8). The inclusion of this structure notably alters the predicted crustal deformation field in regions away from plate boundaries (cf. Supplementary Fig. 5a, b), increasing the seismicity prediction power to 88.89%. The predicted crustal deformation shows a more localized pattern with high deformation rates in the WL and at the edge of CP (Supplementary Fig. 5b). Importantly, the 3D crustal effective viscosity structure significantly reduces the coupling between the North American plate and the fast-moving Pacific plate, as evidenced by the reduced clockwise rotation in the predicted surface motion (Supplementary Fig. 5f). This decoupling effect reduces the average residual to 2.53 mm/yr, while the average angular misfit (43.4°) remains similar to M7.

When further adding basal tractions (M10) to M8, the predicted crustal deformation field notably changes in the interior of WUS, particularly in the SRP and northern RM (cf. Supplementary Fig. 5b, d). The increased crustal deformation rate in the interior of WUS corresponds to a seismicity prediction power of 96.94%. The predicted crustal motion generally exhibits more eastward motion compared to M8 (cf. Supplementary Fig. 5f, h, magenta arrows) due to the eastward drag exerted by the mantle flow, with a similar average residual (2.53 mm/yr) and a reduced average angular misfit (34.0°). When compared to M9, the inclusion of the 3D crustal structure (M10) significantly reduces the magnitude of surface velocity residual while leaving the average angular misfit unchanged.

Effects of lithosphere–asthenosphere interaction on crustal rheology
One notable rheological feature in the geophysically inferred 3D lithospheric effective viscosity structure is a long crustal weak zone along the lithospheric thickness step (Fig. 2g). This weak zone is important in reproducing the fast-deforming ISB (Fig. 5 and Supplementary Fig. 5). Its effective viscosity is independently inferred and is likely controlled by the long-term thermal, compositional, and mechanical evolution on crustal rheology.

Some studies suggested that this weak zone is a result of local volcanism since the Cenozoic due to the intruding hot material6 or lithospheric delamination87. However, a recent geodynamic simulation indicated that there has been prominent hot asthenospheric material intruding from the oceanic mantle beneath the whole WUS since ~12 Ma36. This broad intrusion would have heated up the WUS lithosphere uniformly rather than generating localized weak zones. While the elevated lithospheric GPE resulting from high plateaus may assist the crustal deformation and weakening around the CP6,8, the reconstructed high Nevadaplano2 is much broader than the narrow weak zone. This situation is similar to the current WUS, where lithospheric GPE gradients drive broad crustal deformation7,9 (cf. Fig. 5a, c). Therefore, we propose that the long-term strain localization caused by lithosphere–asthenosphere interaction along the lithospheric thickness step represents a more plausible mechanism for the formation of localized intraplate weak zones and the subsequent development of intraplate seismic belts.

Formation of the lithospheric thickness step
The formation of the lithospheric thickness step likely depends on the contrasting physical properties of the cratonic and non-cratonic lithosphere, as well as the ancient tectonic processes in the WUS, such as lithospheric delamination87 and flat subduction88. Further research involving time-dependent geodynamic simulations that concurrently compute lithospheric and convective mantle dynamics in a unified physical frame is necessary to investigate how the lithospheric thickness step formed and influenced long-term tectonic processes in the WUS.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Summary\nThe manuscript argues that the complicated lithospheric structures of the WUS are the\nprimary cause of the region’s widespread and intense intraplate deformation and uses\nextensive geodynamic simulations to support the argument. The authors first present a\nmodel with realistic lithospheric structures capable of producing crustal deformations that\nmatch the observations significantly better than models with lower degrees of\ncomplexity. They then perform experiments by varying each feature in the model\nindependently to test its effect on the predicted crustal deformations. They finally discuss\nthe possible connections between recent volcanisms in the WUS and the deformations\ncaused by the heterogeneous lithospheric structures in their models.\nRecommendation\nI personally find the results of the paper interesting and are potentially suitable for\npublishing on Nature Communications because they may provide additional insights into\nthe mechanisms of the WUS’s complicated intraplate deformations, an outstanding\nquestion in solid-earth geosciences. Nonetheless, in my opinion, the paper in its current\nform is poorly organized, has major issues with its input parameters, and contains\nnumerous typos and grammatical errors, which call for an overhaul of the manuscript. I\nthus recommend a major revision and am happy to review it again after the resubmission.\nPlease see my detailed reviews below.\nMajor issues\nThe most outstanding issue with the paper is its organization. The contents of the sections\n“Reproducing Lithospheric and Mantle Dynamics” and “Effects of the Heterogeneous\nLithosphere on Crustal Deformation”, two main sections presenting their results, are\nsimilar and their relation unclear. The section “Effects of the Heterogeneous Lithosphere\non Crustal Deformation” lacks a clear focus and contains too many details of too many\nmodels, making it very tedious and difficult to read. I suggest the authors restructure the\ntwo sections to present fewer models in a more organized fashion. In the first section, the\nauthors should show the results of their preferred model, which presumably contains all\nthe realistic heterogeneities in the lithosphere, and compare its crustal-deformation\npredictions with those from a minimalistic model with only the plate boundary\nconstraints and the convective mantle to establish the importance of accounting for\ncomplexities in the lithosphere (The current first section presents three models, whose\nrelations are unclear to me). The second section will then describe the parameter test to\nfind the relative importance of all the complexities. Given this paper is about\ncomplexities in the lithosphere, I feel that the authors only need to vary the viscosity and\ndensity structures of the crust and mantle lithosphere (no need to test the effects of the\nconvective mantle), so they only need to present six models in this section. The authors\nshould also focus on the improvements to the model predictions when complexities are\nadded instead of describing every detail of their results.\nI am skeptical about the lithosphere-thickness model the authors use and the way it is\npresented in the paper. Based on my understanding, the authors use the lithospherethickness may from\nLITHO-1.0, a very coarse-resolution global map estimated mainly\nusing surface-wave data. The authors claim that the map shows significant lithospherethickness\nvariations within the WUS, but the only variation truly within the region is the\none around the Colorado Plateau, whereas the rest of the WUS largely has a constant\nLAB depth of ~65 km, (Fig. 3a). If the authors are to continue using the current LABdepth model, I\nsuggest that they instead describe it as the LAB depth variation between\nthe Colorado Plateau and the rest of the WUS, which is really what it is.\nIn contrast to the coarse-scale LAB-depth map from LITHO-1.0, recent studies using\nbody-wave scattered phases, particularly Liu and Shearer, (2021), have revealed finer\nscale LAB-depth variations within the WUS, which correlate well with the regional\ntectonics. For instance, Liu and Shearer, (2021) showed significantly thinner lithospheres\nbeneath the boundaries of the Basin and Range Province compared with its interior,\nwhich likely will cause additional strain localizations if incorporated into the authors’\nmodels. I understand that running their models using a different LAB model may be too\nmuch of a burden for the authors, but they should at need acknowledge the presence and\npotential effects of finer-scale LAB-depth variations in the WUS instead of ignoring the\nrecent significant progresses in seismically characterizing the LAB. In fact, Fig. S3\nshows that compared to the predictions of their preferred model, the observed\ndeformation rate shows a significantly greater degree of localization along the east and\nwest boundary of the northern Basin and Range Province, the exact locations where Liu\nand Shearer, (2021) observed local thinning of the lithosphere. Therefore, using a finerscale LAB-depth\nmodel will likely further improve the agreement between their model\npredictions and the observations.\nI am also concerned about the authors’ input lithospheric density models. First, the\nmethod section does not clearly state if the density in the crust is spatially variable\nbesides the difference between the Snake River Plain and the rest of the WUS. Second,\nthe authors appear to have modified the thicknesses of the crust and lithosphere based on\nthe surface topography by assuming Airy Isostasy without explicitly stating it. I thus\nsuggest that they clarify if Airy Isostasy is indeed assumed. Third, I am concerned if the\nauthors’ modified crustal-thickness map is significantly different from the observed one,\nwhich is relatively well-constrained by seismological studies. I thus suggest that they\nshow a comparison between their input crustal thickness and the seismologically\nconstrained one in the supplementary material. Lastly, given the importance of the input\ndensity models, I suggest that the authors describe how they built them in greater detail in\nthe method section instead of referring to a previous work unless the journal explicitly\nlimits the section’s length.\nRegarding the discussion on recent volcanisms at the end of the paper, I do not think their\nresults support their claim that a connection exists between the volcanisms and the\nlithosphere extensions predicted by their model. Figs 10b and c show that most of the\nBasin and Range Province is in extension and there is no significant difference between\nits interior and boundaries. Nonetheless, the volcanisms are predominantly located near\nits boundaries, suggesting that their distribution is controlled by processes not captured\nby the model, such as the localized lithosphere thinning at the boundaries of the Basin and\nRange Province shown by Liu and Shearer, (2021). Therefore, I suggest deleting the\ndiscussion on recent volcanisms.\nThe manuscript also contains numerous typos and grammatical errors, some of which are\nlisted below. Although they generally do not cause difficulty in evaluating the merit of\nthe paper, they do leave me the impression that the authors are uncareful while preparing\nthe manuscript.\nMinor issues\nLine 16: The juxtaposition of the terms “lithospheric” and “mantle” is inappropriate as\nthe lithosphere also includes a mantle part. The authors should use “convective mantle”\ninstead of “mantle”. They made this mistake at several other places, e.g., in Line 106.\nLine 31: The authors should use the past tense while describing previous works. For\nexample, here, “realize” should be “realized”. This applies to all other places where\nprevious studies are mentioned.\nLine 68: What is considered “intraplate volcanism” by the authors? Specifically, are the\nCascadia arc volcanisms included? Fig. 1b seems to suggest that they are. Please clarify.\nLine 134: The term “Mesozoic Farallon slab” is confusing. My understanding of the\nmodeling process is that they only compute a present-day snapshot of the mantle-flow\nfield instead of running the model from the Mesozoic. Does it mean the Mesozoic\nFarallon plate that has since been subducted to its current location beneath the WUS? If\nso, I would call it the “subducted Farallon slab” to avoid confusion.\nLines 141–144: How does the model “reproduces the late-Cenozoic volcanic history and\ntopographic changes”? How do they predict the SKS splitting measurements from the\nmantle flow and density anomalies?\nLine 184: Reference 47 is not about seismic attenuation and thus probably a typo.\nLine 215: The sentence should be changed to “…which may have led to an\nunderestimation of …”\nLine 219: “earthquake frequency” should be changed to “seismicity” as the latter is more\ncommonly used in the seismology community. Besides, why does the second invariant of\nstrain rate reflect the deformation rate?\nLines 222–224: In Fig. 4a, why does the area around the Wyoming craton have a very\nlow deformation rate?\nLine 238: The description of M2 does not seem to match that in Table S1 because here,\nM2 is described to have a constant density, whereas in Table S1 it has “the preferred\ndensity structure”. This is related to my aforementioned confusion about whether the\ncrustal density is allowed to vary laterally.\nLine 264: The authors mention changing the “definition of lithospheric GPV” multiple\ntimes in the manuscript, but how can they change the definition of a physical quantity?\nLine 274: The word “linearly” is too strong. Use “monotonically” instead.\nLines 279–281: Can the authors explain why the subduction of the Juan de Fuca plate\ncauses northward motion of the Pacific Northwest in M4, which seems to be roughly\nperpendicular to the subduction direction?\nLine 317–318: I do not see a significantly increased strain localization along the\nlithosphere step in Fig. 5c compared to 5b.\nLine 320: What “to the east of the lithosphere step” means is unclear. I guess it really\nmeans to the east of the boundary between the Basin and Range Province and Colorado\nPlateau. See my previous criticism about the authors’ overstatement about the\nlithosphere-thickness variation in their model.\nLines 330–332: According to Table S1, M7, M8, and M9 have the same structure in the\nconvective mantle but different lithospheric viscosity structures, which is opposite to\nwhat this sentence states.\nLine 336: In which areas does the subducting plate increase crustal deformation?\nLine 443–445: It is impossible to visually determine if the predicted stress field is\nconsistent with the observed focal mechanisms from Fig. 9a. I suggest either removing\nFig. 9 and the associated discussions or making a new figure showing the predicted stress\nfield and the stress field estimated from the focal mechanisms using stress inversion.\nBesides, Fig. 3 does not have Panel c or e, so it’s probably a typo.\nLine 468: Fig. 3 does not have Panel c.\nLines 506–508: This statement is too strong. Fig. 5 and 6 show that the crustal-viscosity\nvariation is at least as important as the lithosphere-thickness variation.\nLine 512–514: In fact, among the seismically active regions, only the ISB is located near\nthe LAB step. Furthermore, even the expected high deformation rate of this region can be\npartially reproduced using models with only the crustal-viscosity heterogeneities (e.g.,\nM8).\nLines 524–526: I don’t think the authors can preclude the local-volcanism origin of the\ncrustal weak zone based on the theory cited here because it represents only one possible\nscenario of hot-material intrusion. Besides, “will heat” should be “would have heated.”\nLines 527–531: I don’t understand the model proposed here.\nLine 585: “allow” should be “allowing”.\nLines 589–592: What are the boundary conditions on the northern and southern\nboundaries of the WUS, which are not bounded by the ocean or the CEUS?\nLine 601: “4” should be a superscript.\nLines 608–610: What are the effects of the MLD on the density and viscosity of the\nmantle lithosphere?\nLines 612–614: This sentence is grammatically incorrect. The authors should either add\nan “and” after the comma or change the comma to a semicolon.\nLine 626: A “the” is missing before “3D”.\nLines 636–637: The sentence starting with “Then” is not a complete sentence.\nLines 653–654: A “the” is missing before “3D”.\nLine 670: A “the” is missing before “vectorial”.\nFig. 3: Please add the legends showing what the two sets of arrows represent in (b).\nBesides, the caption should emphasize that the azimuthal seismic anisotropy in (a) is\nobserved, not predicted, which confused me for quite a while.\nFigs. 4d–f: Add legends.\nFigs. 5d–f: Add legends.\nFig. 9: The green dashed lines in (b) can be mistaken as the green bars denoting\nextensional stresses. Consider change the color.\nPlease add longitude and latitude labels to all maps.\nTable S1: What does “Mantle structure” mean? I assume it means structures in the\nconvective mantle. Besides, “empty” is inappropriate. Use “none” instead.", "reviewer_2": "The paper “Western US Intraplate Deformation Controlled by the Complex Lithospheric Structure” by\nCao & Liu uses numerical modeling to investigate the contributions of various forces on deformation and\nthe distribution of recent volcanism throughout the US. This paper shows somewhat quantitively that an\nunderappreciated control on magmatism, deformation, and earthquake distribution in the western US is\nrelated to gradients in lithospheric structure.\nThe authors should be applauded for taking on such an ambitious project, which brings together\ngeodynamic modeling, and geologic and geophysical observations. I think the major conclusions of the\nauthors as stated above are valid, however their reasoning and clarity in communicating how they\nreached that conclusion is lacking. Much of this is due to the vast amount of models with different\nparameters that they analyze, which is extremely difficult to follow. This could be fixed by adding a table\nof M1-M9 with the different conceptual parameterizations (e.g., columns of crustal thickness, LAB\nthickness, rheology, topography, and rows of “constant” or “variable”). This would make it much easier\nto follow the comparisons. In addition, comparisons between model fits to data are fundamentally\nqualitative in how they are currently presented, making it difficult to know if the conclusion reached by\nthe authors is quantitatively supported by their modeling. Another column could be added to the\nencouraged additional table for misfit to various datasets, like strain, stress orientation, and GPS\nmeasurements, even if they only represent relative misfits compared to each other rather than to true\nvalues, since many of these are only weakly constrained.\nOther issues in clarity exist due to semantic definitions (e.g., L84 “the mantle and lithospheric dynamics”.\nI assume mantle convection was meant here since mantle can also be lithosphere. Another example is\nbroadly using basal traction to apparently refer to any mantle convective force between the\nasthenosphere and lithosphere, rather than purely referring to surface-parallel forces). In a framework of\nunderstanding contributions of difference forces, clearly defining these forces is important. This would\nalso help with the fact that the authors really don’t make any clear conclusions about the relative\ncontributions of body forces, basal tractions, or plate boundary forces, rather simply stating in the\ndiscussion that the lithospheric thickness gradients are important. It isn’t clear which force lithospheric\nthickness gradients should be associated with (it seems that this would affect both body forces and\nviscous drag, or the authors preferred term, “basal tractions”). Thus, the authors never clearly tie their\nintroduction with their results and discussion.\nRegarding the methods, I think there are some unclear points. The authors state that they create a\n3D viscosity model “adjusted to match the observed direction patterns of crustal stress”. This\nseems very circular, as they then go on to say that these viscosity variations are important in\ncontrolling crustal stress orientations, which of course makes sense if they were derived from the\nstress orientations themselves. I think it is important for the authors to clearly identify what is an\ninput to their model vs what is an output. This is exemplified by the L520 “One notable feature along the lithospheric step is the presence of a long weak zone within the crust”. Here, the authors make it sounds like a model result, when it is clearly embedded in the input model.\nIn conclusion, the takeaway of this paper is important, but the strength of their article is not clear in its\ncurrent form due to the lack of quantitative comparisons between the different models (and\ntherefore, plate driving forces, which they discuss quite a bit in the Introduction to set up their\nproblem). Also, it is extremely difficult to keep track of all the models and their different\nassumptions, inputs and results.\nLastly, and perhaps most importantly, the lack of clarity in writing and terminology will prevent this\nmanuscript from having the impact that it likely could have given the importance of the result. My\nrecommendation is that a significant amount of manuscript revision needs to be done to make the paper\nmore intelligible, as it is extremely difficult to read and understand in its current form.\nDetailed Comments:\nAbstract\nL21: “Particularly, the interaction between the mantle flow and lithospheric thickness step along the\neastern boundary of Basin and Range represents a key driving mechanism for localized intracontinental\ndeformation.” Given that the motivation in the Intro is revealing the contributions of plate boundary\nforces, lithospheric body forces, and viscous drag to deformation, which category does this fall under? I\ncould see this being a lithospheric body force (normal force exerted on LAB topography from flow)\nand/or a viscous drag. It’s not clear to me that this study actually illuminates the contributions between\nthe stated body forces. I think the motivation should be changed to clearly apply to the conclusions, or\nthe conclusions should be put in the context of the motivation.\nIntroduction\nL26: uplifting -> uplift\nL26-28: The first part of this sentence is grammatically unclear. I assume you mean that geological and\ngeophysical evidence show clear signs of extension.\nL33-34: This sentence needs to be linked clearly to the stated 3 categories of driving forces. Alternatively,\nyou could walk back the focus on the 3 categories, and simply state that you are exploring something\nthat is largely unexplored in the context of these numerical models.\nL42-44: There has been a lot of research done on the controlling forces of this change in plate motion,\nand is largely attributed to plate boundary forces. An example paper exploring these forces should be\ncited (perhaps a geodesy paper).\nFig. 1: I recommend adding the boundaries of the Northern, Central, and Southern Basin and Range\nprovinces since you discuss them later in the paper. You should also cite the reference for LAB depth\nL68: Need to clarify that PNW volcanoes are not intraplate. The main edifices are typical calc-alkaline\nsubduction-related volcanoes and shouldn't be included with intraplate volcanoes.\nL83-84 and 112: Could you clarify what you mean by \"one-way\" coupling? Seems like a contradiction.\nDoes \"one-way\" coupling mean all forces are drag from asth to lithosphere? What effect does \"two-way\"\ncoupling actually have? It seems that all it implies is that it allows for more spatial complexity in your\nmodel.\nL84, L97", "reviewer_105": "Recommend changing from “lithospheric and mantle dynamics” to “lithospheric and\nmantle convective dynamics” or “crustal and mantle dynamics”\nL130-L138: There don’t seem to be any geological constraints described in this paragraph. All are\ngeophysical\nL140: This language makes it very difficult to understand what is a constraint and what is a result. Do you\nforce your model to fit topographic changes? How do you do this as a function of time? What do you\nmean your model reproduces volcanism?\nL149-150: Isn't the slab from seismic tomography, rather than consistent with it? How does it get\nmodeled as segmented?\nFig 3: It would be very helpful to quantify the fit between SKS and predicted flow if you are going to use\nthis as justification for a \"better\" model. In addition, why should azimuthal anisotropy (an absolute\nmeasure of deformation/flow) correlate with predictions relative to the NA reference frame? I know this\ncomparison is done a lot in geodynamic modeling studies, but it isn’t clear to me why it would be valid,\nas the data and results are in different reference frames\nL211-216: This sentence seems out of place and should probably be in the Introduction/Motivation\nsection. It doesn't seem to particularly contribute to the discussion in this paragraph. Also might help\nwith clarity if it is broken into multiple sentences.\nL226: dominant-dominated\nFig. 4: Very nice conceptual example of the effects of the different parameterizations. I would\nrecommend this difference be quantified somehow to show the point quantitatively. Also, change\ncaption (b) for clarity: “a lithosphere of constant thickness but varying Moho depth” rather than “same as in M1\nL260: This section title (“Effects of the heterogeneous lithosphere on crustal deformation”) is trivially\ndifferent from the prior section ( “Mantle and crustal deformation controlled by the complex lithospheric\nstructure”), neither of which really clearly state the parameters you are going to test in these sections.\nAgain, clearly defining the differences in the models through a table or some other relatively simple\nfigure would greatly help the clarity of which models are testing which contributions.\nL285: May be worth noting briefly how this viscosity is constrained since it seems to have a\nsignificant effect (and/or referencing the methods)\nL311: Which lithospheric layers are varying? Both crustal thickness and total lithospheric\nthickness? Also how does this vary from M3? Clearly the results are different but it isn’t easy to see\nhow/why it varies from M3\nL319: This seems contrary to your results (Fig.5b,c). The model that has variable viscosity but a\nuniform lithosphere (M5) reproduces very similar predicted strain magnitudes to the varying\nlithospheric structure (M6). Perhaps quantifying misfit would make your point more clear, but it\ncertainly isn't qualitatively clear that it makes a big difference in stress and strain predictions\n(although seems to make a difference for GPS predictions)\nL326: Should this be a new subsection? Frankly, it isn’t clear in this paper when we are transitioning from\na Section to a Subsection, and these should be clearly indicated in the section headers\nL338-340: Are you talking about only the black arrows here? Should point that out in the reference to\nthe figures (e.g., cf. 5d, 6d, black arrows)\nL367: I also find the comparison interesting between M7 and M8. We can clearly see that the magnitude\nof surface motion is strongly linked to crustal rheology, however the orientations remain nearly identical.\nHowever, this comes back to the fact that is glossed over in the methods, where your viscosity field is\nderived from surface motions. Regardless, this is another place that would benefit from a quantitative\nanalysis rather than purely qualitative comparison\nL377-380: This statement is unclear to me. It sounds like the significance of the interaction is that it is not\nsignificant given no change in local stress. Maybe change to “variations in the viscosity of the convective\nmantle clearly plays a secondary role relative to LAB gradients in controlling crustal stresses”. Given this,\nit would imply that viscous drag (or at least basal traction) is a negligible plate driving force in the\ninterior. This would be good to link to the overall discussion on driving forces.\nL380-382: May be worth labeling this on the figures? I don't see particularly large changes in defm rate\nbetween the two models.\nL386-389: I assume these results are for M9? I think it would help to explicitly state somewhere which\nmodel these are the predictions for. I assume M9 but it isn't stated. It is hard to keep track of all these\nmodels and their differences\nFigure 7: State which model this is for: M7? This is also a place where an explicit definition of traction\nwould be helpful. I would assume that (a-c) is the traction parallel to the horizontal plane, but there are\nclear upward and downward tractions being observed, which I would think should be called normal\nforces. Again, maybe a semantical difference.\nL435: Change “lithospheric thickness and strength” to “lithospheric properties”. As you modify\ncrustal thickness as well as lithospheric thickness, the former statement doesn’t clearly show all the\nfactors you tested. The “properties” is more vague\nL440: “pressure drag” is this a basal traction? or a body force? This isn’t made clear\nL440: Should this be Fig 7? It isn't stated that lithosphere changes at all in Fig 8.\nL447: Should this be Fig 9a?\nL460-461: This probably should be cited if you are going to state this. Given the long active tectonic\nhistory of the western half of North America, many of these faults are reactivated from prior tectonic\nevents and can rupture even if they are not presently at the optimal orientation relative to the driving\nforces since they will be weaker than unbroken rock.\nL465-468: Many authors attribute the N-S compression to SAF-related tectonic forces as well as\nJdF plate boundary forces. This should probably be added. How can you decouple plate boundary\nforces from basal tractions in strike-slip regions?\nL470-: A couple things here:\nL468-473: How does combining N-S extensional stress with W-E extensional stress lead to W-E\nextension? This seems to overstate the importance of basal traction given that it doesn't seem to\nalter GPE predictions. 2) Where are you talking about when you say southern Basin and Range? The\nsouthern Basin and Range has a relatively specific location (central to southern Arizona and Mexico; see\nJones et al 1992), but there are only 2 or 3 focal mechanisms in this area, only 1 of which seems to be\nN-S oriented. 3) I would bet that including a model of lithospheric properties in Mexico (such as the\nsubducting slab and LAB variations) would significantly alter the results in the southernmost portion of\nthis model (away from the SAF). This should probably be acknowledged and likely makes the discussion\non what factor is controlling stress orientation at these southern latitudes only loosely constrained by\nmeasurements in the US\nL494-496: But which of the 3 forces is this? You set up a discussion that states that this study\ncontributes to answering which dominates the WUS, but never explicitly link the results to the forces\nL508-509: definition of GPE -> do you mean the magnitude of GPE? The definition of GPE doesn't\nchange, it is the values included in the calculation that change.\nL509-512: remove also: it seems that the big change between how the lithosphere and asthenosphere\nact IS through changing the lateral viscosity structure, and thus ease of asthenospheric\nconvention/motion\nL516: I recommend using a phrase besides “transition zone” for clarity, like the \"transition from thin to\nthick lithosphere\". The transition zone has definitions to the geophysical community (mantle transition\nzone) as well as the tectonic community in the western US (transition zone is a small tectonic strip\nbetween southern B&R extension and the Colorado Plateau). Both of these are pretty aseismic and may\njust confuse readers\nL520: This feature is embedded in your model through your parameterization and the depth to the\nMoho and LAB. However, you state it here like it is a result of your modeling rather than an input.\nYou should either delete this sentence (the 2nd sentence is the important one for discusssion anyway) or\nmodify it to clearly differentiate results from inputs\nL521: change \"may reflect\" to \"its rheology is controlled by thermal, compositional, and mechanic\nproperties.\" Clearly it must be controlled by these, as the weak zone in your crust is embedded\nthrough your own parameterization of the model\nL549-554: I do not think that the widespread presence of the low velocity zone necessitates non-vertical\ntransport of melt. In your previous paragraph you conclude that (and I agree) that translithospheric\nextension is likely necessary to get asthenospheric melt to the surface. Given the rather small volume of\nmost (non-Yellowstone related) magmatism in the WUS, why couldn’t these localized extension region\nonly tap the asthenosphere below them, while melt away from these zones remains trapped below an\nimpermeable lithosphere due to no vertical strain?\nIn addition, I don't think it is appropriate to compare mostly small-volume basaltic (i.e., very low\nviscosity) magmatism that spatially dominates the intraplate WUS to a long-lived rhyolitic (i.e., high\nviscosity) magmatism, and perhaps even less appropriate to compare locations of short-lived\nasthenosphericly-sourced magmatism to long-lived magma residence in magma reservoirs. I agree that\nbelow more silicic magmatic centers, there is no need for the reservoir to lie below the edifice, as the\nedifice location is likely controlled by brittle, active structures in the upper crust (especially for high\nviscosity magmas). There is plenty of literature on that. But I do not think this applies to pathways for\nvery low viscosity melts. Also, you are implying that melt is able to be taps from 100s of km away from\nwhere they erupt with this comparison, whereas active magmatic systems generally show reservoirs\nbelow edifices, while offset, tend to be no more than ~10 km away from the edifice (Lerner et al., 2020).\nThis mix increase with depth, but it would be quite a stretch to say that the center of the B&R is sourcing\nmagmatism near the Wasatch. Lastly, a simple calculation comparing the volume of the LVZ below the\nB&R with some assumed melt percentage with the total erupted volume in the B&R will show that very\nlittle assumed melt reaches the surface.\nYou state a more plausible interpretation in L554-556, by saying melt reaches the surface if lithospheric\nstresses allow for it. In a more compressional state, the lithosphere could probably be inferred to act as\nan impermeable cap, even if melt is present.\nThis could all be corrected by simply deleting the sentence encompassed by L553-554\nL560-561: It would be nice to have more specific examples rather than just pointing out large regions. I\nthink the East African Rift system near cratons and the Anatolian system are both areas where literature\nhas interpreted that magmatism correlates with throughgoing lithospheric structures and presumed\nchanges in LAB depth\nMethods:\nL607: This is a really low resolution model of crustal thickness. Why not use better constrained\nthicknesses, like from Schmandt et al 2015 (model on Earthscope EMC page)\nL627-629: This seems circular...Of course adding 3D viscosity variations to your model in this study leads\nto better fits to crustal stresses, the viscosity model is derived from fitting the stresses!"}, "author_response": "Below are our detailed point-to-point responses to the reviewers’ comments and suggestions, with\noriginal comments in black and our response in blue.\nTo reviewer #1:\nI personally find the results of the paper interesting and are potentially suitable for\npublishing on Nature Communications because they may provide additional insights into\nthe mechanisms of the WUS’s complicated intraplate deformations, an outstanding\nquestion in solid-earth geosciences. \nNonetheless, in my opinion, the paper in its current form is poorly organized, has major issues\nwith its input parameters, and contains numerous typos and grammatical errors, which call for an\noverhaul of the manuscript. I thus recommend a major revision and am happy to review it again\nafter the resubmission. Please see my detailed reviews below.\nResponse: Thanks for your kind words and insightful comments regarding our manuscript. We\ngreatly appreciate your constructive feedback. To enhance the quality of our manuscript, we have\nmade significant revisions based on your comments and suggestions. Specifically, we have\nreorganized our “Results and Discussion” section to enhance the logical flow and clarity.\nAdditionally, we have expanded the “Methods” section to provide more details about our model\ninputs, assumptions, and justifications for our modeling approach. Furthermore, we have\nintroduced a discussion on the potential links between our model and recent advancements in\nseismological studies and deleted all volcanism-related content. Thanks again for your\nconstructive reviews.\nThe most outstanding issue with the paper is its organization. The contents of the sections\n“Reproducing Lithospheric and Mantle Dynamics” and “Effects of the Heterogeneous\nLithosphere on Crustal Deformation”, two main sections presenting their results, are\nsimilar and their relation unclear. The section “Effects of the Heterogeneous Lithosphere\non Crustal Deformation” lacks a clear focus and contains too many details of too many\nmodels, making it very tedious and difficult to read. I suggest the authors restructure the\ntwo sections to present fewer models in a more organized fashion. In the first section, the\nauthors should show the results of their preferred model, which presumably contains all\nthe realistic heterogeneities in the lithosphere, and compare its crustal-deformation\npredictions with those from a minimalistic model with only the plate boundary\nconstraints and the convective mantle to establish the importance of accounting for\ncomplexities in the lithosphere (The current first section presents three models, whose\nrelations are unclear to me). The second section will then describe the parameter test to\nfind the relative importance of all the complexities. Given this paper is about\ncomplexities in the lithosphere, I feel that the authors only need to vary the viscosity and\ndensity structures of the crust and mantle lithosphere (no need to test the effects of the\nconvective mantle), so they only need to present six models in this section. The authors\nshould also focus on the improvements to the model predictions when complexities are\nadded instead of describing every detail of their results.\nResponse: Thanks for your valuable suggestions! We have divided the “Results and Discussion”\nsection into several subsections to enhance the logic flow. In the revised “Results and Discussion”\nsection, we first present how we construct the fully dynamic 3D model. Then, we discuss the\nimpact of lateral varying lithospheric thickness on the asthenospheric flow. Following that, we\nshow our results on crustal deformation in two sections. In the first section, we show the necessity\nof a complex 3D lithospheric structure in reproducing the observed crustal deformation pattern.\nThen, we show how a heterogeneous lithosphere influences the contribution of mantle convection\nto crustal deformation through the plate boundary effect and basal tractions and estimate the crustal\ndeformation and stress resulting from lithospheric body forces. After that, we present the results\non locally increased lithosphere-asthenosphere interaction beneath the western US (WUS). Finally,\nwe discuss the relative importance of the proposed driving mechanisms accounting for the realistic\nlithospheric structure.\nPrevious studies using simplified lithospheric structures suggested that the mantle convection (e.g.,\nSteinberger et al., 2001, EPSL; Becker et al., 2015, Nature) and crustal gravitational potential\nenergy (e.g., Jones et al., 1996, Nature; Artyushkov, 1973, JGR) are two main drivers for crustal\ndeformation. Therefore, we consider these two factors first in the “Complex Lithospheric Structure\nControls WUS Crustal Deformation” section. However, the predicted crustal deformation field\nfails to reproduce the observed crustal deformation field. Then, we show the predicted crustal\ndeformation from our best-fitting model to highlight the importance of a complex three-\ndimensional lithosphere.\nIn the “WUS Crustal Deformation due to Different Driving Mechanisms” section, we have\nretained all six models because lateral variations in the lithospheric thickness have a profound\nimpact on the interaction between the continental lithosphere and convecting mantle, including the\nplate boundary effect and basal tractions. A heterogeneous lithosphere will impact how the plate\nboundary forces transmit into the interior. Also, when there is a flat lithosphere-asthenosphere\nboundary (LAB), the convecting mantle exerts horizontal shear and normal force in the radial\ndirection along the LAB. However, in regions with a strong lateral gradient in lithospheric\nthickness, such as along the eastern boundary of the Basin and Range (B&R), the mantle flow\nexerts normal force in the horizontal direction, as demonstrated in Fig. 7e and 7f. This force was\npreviously overlooked in the studies assuming a flat LAB (e.g., Flesch et al., 2000, Science; Becker\net al., 2016, Nature; Bahadori et al., 2022, Nature Communications). Therefore, we have retained\nall six models to fully explore the importance of a heterogeneous lithosphere on the contribution\nof mantle convection to crustal deformation. To enhance the clarity in this section, we have revised\nthe text to focus on the predicted crustal deformation patterns and the improvements when\ncomplexities are added, according to your suggestions.\nI am skeptical about the lithosphere-thickness model the authors use and the way it is\npresented in the paper. Based on my understanding, the authors use the lithosphere-thickness may from\nLITHO-1.0, a very coarse-resolution global map estimated mainly\nusing surface-wave data. The authors claim that the map shows significant lithosphere-thickness\nvariations within the WUS, but the only variation truly within the region is the\none around the Colorado Plateau, whereas the rest of the WUS largely has a constant\nLAB depth of ~65 km, (Fig. 3a). If the authors are to continue using the current LAB-depth model,\nI suggest that they instead describe it as the LAB depth variation between\nthe Colorado Plateau and the rest of the WUS, which is really what it is.\nResponse: We really appreciate your insight! Our lithospheric thickness model is a hybrid one,\nincorporating LITHO 1.0 (Payanos et al., 2014, JGR) and a regional body-wave tomography model\nby Schmandt & Lin (2014, GRL). LITHO 1.0 characterizes a thin lithosphere with smooth lateral\nvariations in the WUS, including the Intermountain West (IW) and Colorado Plateau (CP). In\ncontrast, the body-wave tomography revealed a continuous fast anomaly extending beyond 150\nkm in the IW and CP. To capture the first-order lithospheric thickness variation, we combined\nLITHO 1.0 with the body-wave tomography model to create a continental-scale lithospheric\nthickness model. In practice, we adjusted the WUS lithospheric thickness derived from LITHO\n1.0 by tracking continuous fast anomalies in the body-wave tomography model. We have added a\nsection in the revised “Methods” describing how we construct the lithospheric thickness model\n(lines 1738-1745). Also, in the new section, “Construction of the Fully Dynamic 3D Model”, we\nhave specified that the lithospheric thickness variation discussed in the manuscript is the\nlithospheric thickness step along the eastern boundary of B&R (lines 444-447).\nIn contrast to the coarse-scale LAB-depth map from LITHO-1.0, recent studies using body-wave\nscattered phases, particularly Liu and Shearer, (2021), have revealed finer scale LAB-depth\nvariations within the WUS, which correlate well with the regional tectonics. For instance, Liu and\nShearer, (2021) showed significantly thinner lithospheres beneath the boundaries of the Basin and\nRange Province compared with its interior, which likely will cause additional strain localizations\nif incorporated into the authors’ models. I understand that running their models using a different\nLAB model may be too much of a burden for the authors, but they should at need acknowledge\nthe presence and potential effects of finer-scale LAB-depth variations in the WUS instead of\nignoring the recent significant progresses in seismically characterizing the LAB. In fact, Fig. S3\nshows that compared to the predictions of their preferred model, the observed\ndeformation rate shows a significantly greater degree of localization along the east and\nwest boundary of the northern Basin and Range Province, the exact locations where Liu\nand Shearer, (2021) observed local thinning of the lithosphere. Therefore, using a finer-scale LAB-depth\nmodel will likely further improve the agreement between their model\npredictions and the observations.\nResponse: Thanks for directing us to recent seismological studies. Using a finer-scale lithospheric\nthickness model may further improve the agreement between our model prediction and the\nobservations. However, we want to highlight the first-order importance of lithospheric thickness\nvariation on lithospheric deformation and shallow mantle flow, as they were neglected in most\nprevious geodynamic studies. Thus, a regional-scale lithospheric thickness model, including the\nfirst-order lithospheric thickness contrast between the actively deforming WUS and the rest of the\nconterminous US, is enough for this purpose. To acknowledge the recent advances in\nseismological studies, we have introduced a paragraph discussing our model results and the\nobserved lithospheric thinning in the revised manuscript (lines 1254-1262).\nI am also concerned about the authors’ input lithospheric density models. First, the\nmethod section does not clearly state if the density in the crust is spatially variable\nbesides the difference between the Snake River Plain and the rest of the WUS. Second,\nthe authors appear to have modified the thicknesses of the crust and lithosphere based on\nthe surface topography by assuming Airy Isostasy without explicitly stating it. I thus\nsuggest that they clarify if Airy Isostasy is indeed assumed. Third, I am concerned if the\nauthors’ modified crustal-thickness map is significantly different from the observed one,\nwhich is relatively well-constrained by seismological studies. I thus suggest that they\nshow a comparison between their input crustal thickness and the seismologically\nconstrained one in the supplementary material. Lastly, given the importance of the input\ndensity models, I suggest that the authors describe how they built them in greater detail in\nthe method section instead of referring to a previous work unless the journal explicitly\nlimits the section’s length.\nResponse: Thanks for the comments! We have added a detailed description of our lithospheric\ndensity structure in the revised “Methods” section (lines 1774-1787). In our preferred lithospheric\ndensity structure, we do not allow lateral density variations in the crust, except for the Snake River\nPlain (SRP) due to its enrichment in basaltic composition (DeNosaquo et al., 2009, JVGR;\nMcCurry & Rodgers, 2009, JVGR).\nDuring the modification of lithospheric and crustal thickness, we did not assume Airy Isostacy.\nPractically, we adjusted the LAB and Moho in two steps, and we ran a simulation with all density\nanomalies in the continental lithosphere and convective mantle in each step to get the predicted\ntopography, which contains both dynamic and isostatic components. Based on the difference\nbetween predicted topography and observed topography, we first adjusted the LAB depth due to\nits larger uncertainty compared to the Moho depth. Then, we did another iteration to adjust the\nMoho depth. The resulting Moho depth is consistent with recent seismological studies. We have\nalso included a figure in the supplementary material (Supplementary Fig. 5) showing the\ncomparison between our input Moho model and recent seismologically imaged Moho in the WUS\n(Laske et al., 2013; Schmandt et al., 2015, GRL).\nRegarding the discussion on recent volcanisms at the end of the paper, I do not think their\nresults support their claim that a connection exists between the volcanisms and the\nlithosphere extensions predicted by their model. Figs 10b and c show that most of the\nBasin and Range Province is in extension and there is no significant difference between\nits interior and boundaries. Nonetheless, the volcanisms are predominantly located near\nits boundaries, suggesting that their distribution is controlled by processes not captured\nby the model, such as the localized lithosphere thinning at the boundaries of the Basin and\nRange Province shown by Liu and Shearer, (2021). Therefore, I suggest deleting the\ndiscussion on recent volcanisms.\nResponse: Thank you for this comment, which we agree with. In the revised manuscript, we have\ndeleted the discussion on recent intraplate volcanism according to your suggestion.\nThe manuscript also contains numerous typos and grammatical errors, some of which are\nlisted below. Although they generally do not cause difficulty in evaluating the merit of\nthe paper, they do leave me the impression that the authors are uncareful while preparing\nthe manuscript.\nResponse: We really appreciate this comment. During the revision, we carefully went over the\ngrammar and checked all spells, hoping the revision is grammatically correct and logically sound.\nMinor issues\nLine 16: The juxtaposition of the terms “lithospheric” and “mantle” is inappropriate as\nthe lithosphere also includes a mantle part. The authors should use “convective mantle”\ninstead of “mantle”. They made this mistake at several other places, e.g., in Line 106.\nResponse: Thank you for pointing this out! We have changed the expression to “lithosphere and\nconvective mantle”. In addition, we have applied this change to all related places.\nLine 31: The authors should use the past tense while describing previous works. For\nexample, here, “realize” should be “realized”. This applies to all other places where\nprevious studies are mentioned.\nResponse: We really appreciate this suggestion. During the revision, we have changed all places\ndescribing previous studies to the past tense.\nLine 68: What is considered “intraplate volcanism” by the authors? Specifically, are the\nCascadia arc volcanisms included? Fig. 1b seems to suggest that they are. Please clarify.\nResponse: Here, we want to use the term “intraplate volcanism” to describe the volcanoes away\nfrom active subduction zones. The volcanoes in the Cascadia arc should not be included in this\ncase. In the revision, we have removed all volcanism-related content from the manuscript.\nLine 134: The term “Mesozoic Farallon slab” is confusing. My understanding of the\nmodeling process is that they only compute a present-day snapshot of the mantle-flow\nfield instead of running the model from the Mesozoic. Does it mean the Mesozoic\nFarallon plate that has since been subducted to its current location beneath the WUS? If\nso, I would call it the “subducted Farallon slab” to avoid confusion.\nResponse: This is a good point. Our present-day convective mantle structure is taken from a time-\ndependent mantle flow model starting from the Late Cenozoic (20 Ma, Zhou & Liu, 2017, G-\ncubed). The term “Mesozoic Farallon slab” means the Farallon Plate has been subducted beneath\nNorth America since the Mesozoic. As recommended, we have changed it to “subducted Farallon\nslab” to avoid ambiguity (line 412).\nLines 141–144: How does the model “reproduces the late-Cenozoic volcanic history and\ntopographic changes”? How do they predict the SKS splitting measurements from the\nmantle flow and density anomalies?\nResponse: These are good questions! As stated above, our convective mantle structure is taken\nfrom a time-dependent mantle flow model starting from the Late Cenozoic (20 Ma, Zhou & Liu,\n2017, G-cubed). The predicted time-dependent mantle flow is used to investigate the spatial\ncorrelation between the intruding hot asthenospheric material and intraplate volcanism (Zhou et\nal., 2018, Nature Geoscience), calculate flow-induced asthenospheric anisotropy and SKS splitting\n(Zhou et al., 2018, EPSL), and quantify the change in dynamic topography (Zhou & Liu, 2019,\nEPSL). These studies successfully reproduced the late Cenozoic WUS volcanic history, present-\nday observed SKS splitting measurements in the WUS, and WUS topographic changes since the\nMiocene. These calculations are outside of this manuscript’s scope, so we invite the reviewer and\nreaders to check our previous publications for more details. We have revised this part to clearly\nstate that our present-day convecting mantle structure is adopted from a time-dependent mantle\nflow model, and the time-dependent convecting mantle structure matches the mentioned\ngeophysical and geological observations, as shown in our previous studies (lines 407-440).\nLine 184: Reference 47 is not about seismic attenuation and thus probably a typo.\nResponse: Yes, this is a typo. The right reference should be 48 in the original order (50 in the\nrevised manuscript), we have corrected this typo.\nLine 215: The sentence should be changed to “...which may have led to an\nunderestimation of ...”\nResponse: Thank you for this suggestion. We have substantially reorganized and revised this\nsection, and this sentence has been removed from the revised manuscript.\nLine 219: “earthquake frequency” should be changed to “seismicity” as the latter is more\ncommonly used in the seismology community. Besides, why does the second invariant of\nstrain rate reflect the deformation rate?\nResponse: We have changed the term “earthquake frequency” to “seismicity” in this sentence and\nall related parts accordingly.\nGenerally, the second invariant of a tensor is a scalar measurement of its magnitude. Specifically\nfor the strain rate tensor, its second invariant measures the total magnitude of the deformation rate,\nincluding dilatation (volume change) and shear. This is a widely used scalar quantity of\ndeformation rate in the geodetic community (e.g., Kreemer et al., 2003, G-cubed), and it has been\nshown to be a good proxy for seismicity (e.g., Bird & Kreemer, 2015, BSSA; Kreemer et al., 2022,\nSRL).\nLines 222–224: In Fig. 4a, why does the area around the Wyoming craton have a very\nlow deformation rate?\nResponse: This is a good observation! In this model, the crustal deformation pattern is dominated\nby the plate boundary effect, which monotonically decays with distance. The Wyoming plateau\n(WP) is situated far from the plate boundaries, so the deformation rate is low within the WC and\nits surrounding area.\nLine 238: The description of M2 does not seem to match that in Table S1 because here,\nM2 is described to have a constant density, whereas in Table S1 it has “the preferred\ndensity structure”. This is related to my aforementioned confusion about whether the\ncrustal density is allowed to vary laterally.\nResponse: M2 has the preferred crustal density structure, in which the crustal density is constant\n(2850 kg/m3), except for the SRP. Since the crust in the SRP is enriched in basaltic composition,\nwe assume the average crustal density in the SRP is denser (2950 kg/m3) than in other places. This\nassumption is consistent with geological and geophysical studies (e.g., DeNosaquo et al., 2009,\nJVGR; McCurry & Rodgers, 2009, JVGR; Zhou & Liu, 2019, EPSL). We have modified the\ndescription here (lines 757-760) and added details about our crustal density structure in the revised\nMethods session to avoid ambiguity (lines 1774-1787).\nLine 264: The authors mention changing the “definition of lithospheric GPV” multiple\ntimes in the manuscript, but how can they change the definition of a physical quantity?\nResponse: This is a great point! Traditionally, the lithospheric gravitational potential energy (GPE)\nis calculated under the assumption of the existence of a global compensation depth such as at 100\nkm (e.g., Ghosh et al., 2012; Becker et al., 2016, Nature; Bahadori et al., 2019, Geosphere).\nHowever, the assumed global compensation depth is questionable and may vary geographically\n(Schubert & Turcot, 2002). Here, we mean the traditional calculation with this implicit assumption\nby using the term “definition of lithospheric GPE”. However, the word ‘definition’ is indeed a bit\nconfusing. In the revised manuscript, we have changed this term to ‘calculation’, and added a\nsection to describe how we calculate three types of driving forces and discuss the differences\nbetween our calculations and those in previous studies (lines 803-835). Particularly, we pointed\nout the problem embedded in the conventional calculation of lithospheric GPE (lines 827-829).\nLine 274: The word “linearly” is too strong. Use “monotonically” instead.\nResponse: We have changed the word “linearly” to “monotonically”. Thanks for the suggestion.\nLines 279–281: Can the authors explain why the subduction of the Juan de Fuca plate\ncauses northward motion of the Pacific Northwest in M4, which seems to be roughly\nperpendicular to the subduction direction?\nResponse: This is a good observation! In M4, the general northwestward to northward surface\nmotion is dominated by the drag from the fast-moving Pacific plate due to the strong coupling.\nThe subducting Juan de Fuca (JdF) plate causes northeastward surface motion in the subduction\ndirection. The combination of the northwestward motion, caused by the right-lateral shear along\nthe Pacific-North America plate, and the northeastward motion, caused by the subducting JdF plate,\nresults in an apparent northward surface motion in the Pacific Northwest.\nLine 317–318: I do not see a significantly increased strain localization along the\nlithosphere step in Fig. 5c compared to 5b.\nResponse: Here, we want to emphasize that the deformation is concentrated along the lithospheric\nthickness step on the thin and weak side (i.e., the eastern boundary of B&R) in M6 (Fig. 5c).\nIn comparison, the crustal deformation diffuses into the interior of CP and WP in M5 (Fig. 5b).\nWe have modified this part to make our point more straightforward by deleting the statement about\nmore localization. Besides, we have added a description clearly pointing to the change in the\ninterior of CP and WP (lines 978-980).\nLine 320: What “to the east of the lithosphere step” means is unclear. I guess it really\nmeans to the east of the boundary between the Basin and Range Province and Colorado\nPlateau. See my previous criticism about the authors’ overstatement about the\nlithosphere-thickness variation in their model.\nResponse: Here, “to the east of the lithospheric step” means to the east of the lithospheric thickness\nstep along the eastern boundary of B&R (the orange dashed line in Fig. 5c). We have added a\ndescription in lines 444-447 to specify that the lithospheric thickness variations discussed in this\nmanuscript is the sharp lithospheric thickness step along the eastern boundary of B&R.\nLines 330–332: According to Table S1, M7, M8, and M9 have the same structure in the\nconvective mantle but different lithospheric viscosity structures, which is opposite to\nwhat this sentence states.\nResponse: This is a good point. Here, we extracted basal tractions by calculating the difference\nbetween two predicted stress fields with the same lithospheric structure but different convective\nmantle structures. For example, the stress field shown in Fig. 6c (now Fig. 5f) is the difference\nbetween M9 and M6, which have the same lithospheric structure but different convective mantle\nstructures (M9 has all density anomalies in the convective mantle, but M6 has an empty mantle).\nWe have rewritten this sentence to improve its clarity (lines 866-868) and revised the figure to\nonly show the estimated crustal stress field due to different mechanisms when the realistic 3D\nlithosphere is included in the model.\nLine 336: In which areas does the subducting plate increase crustal deformation?\nResponse: In the PNW back arc region, i.e., near the western boundary of the SRP. We have added\na more specific description in this sentence (lines 899-901).\nLine 443–445: It is impossible to visually determine if the predicted stress field is\nconsistent with the observed focal mechanisms from Fig. 9a. I suggest either removing\nFig. 9 and the associated discussions or making a new figure showing the predicted stress\nfield and the stress field estimated from the focal mechanisms using stress inversion.\nBesides, Fig. 3 does not have Panel c or e, so it’s probably a typo.\nResponse: Thanks for the suggestion! We have revised Fig. 9 to include two panels comparing our\nmodel predictions with recent geophysical estimations of S direction and Aφ parameter (Lund\nHmax\nSnee & Zoback, 2020, Nature Communications). Here, Fig. 3 should be Fig. 4. We have corrected\nthis typo.\nLine 468: Fig. 3 does not have Panel c.\nResponse: Same as the above, it should be Fig. 4. We have corrected it.\nLines 506–508: This statement is too strong. Fig. 5 and 6 show that the crustal-viscosity\nvariation is at least as important as the lithosphere-thickness variation.\nResponse: We really appreciate this comment! The 3D crustal effective viscosity is important in\nreproducing the observed crustal deformation pattern. However, we analyzed the potential\nmechanisms for creating the localized crustal weak zone and found that the long-term localized\nlithosphere-asthenosphere interaction may ultimately control it. Therefore, we think this statement\nis generally valid. We have brought up the discussion on formation mechanisms for crustal weak\nzones and merged it with our results showing localized lithosphere-asthenosphere interaction into\na new section, “Mechanism for Localized Intraplate Deformation” (lines 1187-1262).\nLine 512–514: In fact, among the seismically active regions, only the ISB is located near\nthe LAB step. Furthermore, even the expected high deformation rate of this region can be\npartially reproduced using models with only the crustal-viscosity heterogeneities (e.g.,\nM8).\nResponse: Here, we refer to “intraplate earthquakes” as the seismicity in the ISB. The crustal\nrheology likely reflects the long-term thermal, compositional, and mechanical effects of tectonic\nevents. By analyzing the potential driving mechanisms, we propose that the crustal weak zone in\nthe ISB was created by the long-term strain localization caused by lithosphere-asthenosphere\ninteraction along the lithospheric thickness step (lines 1222-1252). Therefore, the varying\nasthenospheric flow ultimately controls the seismicity in the ISB. In the revised manuscript, this\nsentence has been removed.\nLines 524–526: I don’t think the authors can preclude the local-volcanism origin of the\ncrustal weak zone based on the theory cited here because it represents only one possible\nscenario of hot-material intrusion. Besides, “will heat” should be “would have heated.”\nResponse: This is a good point. First, the hot material intrusion predicted by Zhou et al. (2018,\nNature Geoscience) is consistent with the surface volcanism trend since 16 Ma. This broad\nintrusion has been covering the entire upper-most asthenosphere beneath the WUS since ~12 Ma\n(Fig. 4 in Zhou et al., 2018, Nature Geoscience), aligning with the onset of volcanism along the\neastern boundary of B&R. Besides, the same time-dependent mantle flow predicted flow-induced\nseismic anisotropy (Zhou et al., 2018, EPSL) and dynamic topography change since the Miocene\n(Zhou & Liu, 2019, EPSL), consistent with observations over the entire WUS. Therefore, this\nscenario seems to be the most plausible one.\nIn addition, the localized crustal weak zone extends from the southern B&R all the way to the\nnorthern Rocky Mountain (RM), but the Cenozoic volcanism only clustered around the edge of\nCP and in the SRP. Notably, the crust has low effective viscosity in the Wasatch Fault Zone and\nnorthern RM, where no Cenozoic intraplate volcanoes exist. The disagreement between the spatial\npattern of the crustal weak zone and intraplate volcanism indicates the volcanic origin is not valid.\nIn contrast, since ~12 Ma, the eastward hot asthenospheric flow has been continuously interacting\nwith the lithospheric thickness step (Zhou et al., 2018, Nature Geoscience), which spatially\ncoincides with the entire crustal weak zone, providing a more plausible mechanism for the\nformation of the weak zone.\nMoreover, the majority of intraplate volcanoes along the eastern boundary of B&R were fueled by\nsmall-volume basaltic melt with low viscosity, which has limited ability to channel through the\nentire lithosphere. Forming these basaltic volcanoes may further require pre-existing lithospheric\nweak zones, rather than creating the weak zones.\nLines 527–531: I don’t understand the model proposed here.\nResponse: Here, we assess another potential formation mechanism for the narrow crustal weak\nzone. The high plateau, such as the Nevadaplano, will locally increase the extensional stress due\nto the elevated crustal GPE (e.g., Jones et al., 1996, Nature). The increased extensional stress may\ncause intense localized crustal deformation, leading to the formation of a crustal weak zone.\nHowever, the reconstructed Nevadaplano is much wider than the weak zone converted from\nseismic attenuation, and it only covered the boundary between the B&R and CP. Thus, it is hard\nto explain the formation of the long, narrow crustal weak zone extending all the way to the northern\nRM by this hypothesis.\nLine 585: “allow” should be “allowing”.\nResponse: We have changed the word accordingly (line 1600).\nLines 589–592: What are the boundary conditions on the northern and southern\nboundaries of the WUS, which are not bounded by the ocean or the CEUS?\nResponse: This is a good question! Our model includes the conterminous US and parts of Canada\nand Mexico. In the model, we apply plate motions in the surrounding ocean basins, stable central\nand eastern US (east of 110˚E) and Canada (north of 55˚N). Other places all have a free-slip\nboundary condition. We have added some details about the hybrid surface velocity boundary\nconditions in the revised “Methods” section (lines 1605-1608).\nLine 601: “4” should be a superscript.\nResponse: Yes, it is a superscript, and we have corrected it.\nLines 608–610: What are the effects of the MLD on the density and viscosity of the\nmantle lithosphere?\nResponse: We take the MLD as a boundary separating the layered continental lithospheric mantle\n(CLM). In our lithospheric density model, the upper and lower CLM have different compositional\ndensities. The lower CLM is compositionally denser than the upper CLM. This configuration is\nsupported by our recent works (e.g., Zhou et al., 2019, EPSL; Hu et al., 2019, Nature Geoscience;\nCao & Liu, 2021, JGR; Wang et al., 2021, GRL; Wang et al., 2021, JGR; Wang et al., 2023, Nature\nGeoscience). The MLD does not affect the effective viscosity in our model. We have added details\nregarding the effects of the MLD in lines 1744-1772.\nLines 612–614: This sentence is grammatically incorrect. The authors should either add\nan “and” after the comma or change the comma to a semicolon.\nResponse: We have added an “and” after the comma.\nLine 626: A “the” is missing before “3D”.\nResponse: We have added the article before “3D”.\nLines 636–637: The sentence starting with “Then” is not a complete sentence.\nResponse: We completed this sentence by adding “we” before “determine”.\nLines 653–654: A “the” is missing before “3D”.\nResponse: We have added the article in the sentence.\nLine 670: A “the” is missing before “vectorial”.\nResponse: We have added the article before “vectorial”.\nFig. 3: Please add the legends showing what the two sets of arrows represent in (b).\nBesides, the caption should emphasize that the azimuthal seismic anisotropy in (a) is\nobserved, not predicted, which confused me for quite a while.\nResponse: Thank you for the comment! We have modified the figure to include the suggestions.\nFigs. 4d–f: Add legends. Figs. 5d–f: Add legends.\nResponse: In these figures, we include the velocity scale at the bottom of the figure and specify\nthe meaning of the colors in the figure caption.\nFig. 9: The green dashed lines in (b) can be mistaken as the green bars denoting\nextensional stresses. Consider change the color.\nResponse: We have changed the color of the dashed line to orange.\nPlease add longitude and latitude labels to all maps.\nResponse: We have added the longitude and latitude labels to all maps.\nTable S1: What does “Mantle structure” mean? I assume it means structures in the\nconvective mantle. Besides, “empty” is inappropriate. Use “none” instead.\nResponse: Yes, it should be the convective mantle structure. We have revised the supplementary\ntable S1 to address this issue and changed “none” to empty in the table.\nTo reviewer #2:\nThe authors should be applauded for taking on such an ambitious project, which brings together\ngeodynamic modeling, and geologic and geophysical observations. I think the major conclusions\nof the authors as stated above are valid, however their reasoning and clarity in communicating how\nthey reached that conclusion is lacking. Much of this is due to the vast amount of models with\ndifferent parameters that they analyze, which is extremely difficult to follow. This could be fixed\nby adding a table of M1-M9 with the different conceptual parameterizations (e.g., columns of\ncrustal thickness, LAB thickness, rheology, topography, and rows of “constant” or “variable”).\nThis would make it much easier to follow the comparisons. In addition, comparisons between\nmodel fits to data are fundamentally qualitative in how they are currently presented, making it\ndifficult to know if the conclusion reached by the authors is quantitatively supported by their\nmodeling. Another column could be added to the encouraged additional table for misfit to various\ndatasets, like strain, stress orientation, and GPS measurements, even if they only represent relative\nmisfits compared to each other rather than to true values, since many of these are only weakly\nconstrained.\nResponse: Thank you so much for your kind words regarding our work! The comments really\nencouraged us. We have modified the table in the supplementary materials (supplementary table\nS1) to include details of our model setups and quantitative comparisons between our model\nprediction and observation. We hope this revised table with quantitative evaluations could make\nour conclusions easier to be assessed. In addition, we have included several sections in the revised\n“Methods” section to describe how we did the quantitative evaluations between the model\nprediction and observation.\nOther issues in clarity exist due to semantic definitions (e.g., L84 “the mantle and lithospheric\ndynamics”. I assume mantle convection was meant here since mantle can also be lithosphere.\nAnother example is broadly using basal traction to apparently refer to any mantle convective force\nbetween the asthenosphere and lithosphere, rather than purely referring to surface-parallel forces).\nIn a framework of understanding contributions of difference forces, clearly defining these forces\nis important. This would also help with the fact that the authors really don’t make any clear\nconclusions about the relative contributions of body forces, basal tractions, or plate boundary\nforces, rather simply stating in the discussion that the lithospheric thickness gradients are important.\nIt isn’t clear which force lithospheric thickness gradients should be associated with (it seems that\nthis would affect both body forces and viscous drag, or the authors preferred term, “basal\ntractions”). Thus, the authors never clearly tie their introduction with their results and discussion.\nResponse: Thanks for the comments. We agree that the terminology in our original manuscript\nwas loosely defined, and how the lithospheric thickness variations contribute to the driving forces\nwas not clear. In the revised manuscript, we have added a paraph describing our definitions of\ndifferent forces when the realistic 3D lithospheric structure is considered (lines 803-835). In the\nsame paragraph, we also discussed how lithospheric thickness variations impact the proposed\ndriving forces. In addition, we have substantially revised our manuscript to strengthen the link\nbetween our “Introduction” and “Results and Discussion” sections. Particularly, we have added a\nquantitative evaluation of the relative importance of three proposed driving forces in driving the\nWUS crustal deformation (Fig. 9).\nRegarding the methods, I think there are some unclear points. The authors state that they create a\n3D viscosity model “adjusted to match the observed direction patterns of crustal stress”. This\nseems very circular, as they then go on to say that these viscosity variations are important in\ncontrolling crustal stress orientations, which of course makes sense if they were derived from the\nstress orientations themselves. I think it is important for the authors to clearly identify what is an\ninput to their model vs what is an output. This is exemplified by the L520 “One notable feature\nalong the lithospheric step is the presence of a long weak zone within the crust”. Here, the authors\nmake it sounds like a model result, when it is clearly embedded in the input model.\nResponse: Thank you for the constructive comments! Our preferred lithospheric effective viscosity\nstructure was converted from a recent seismic attenuation study (Hearn, 2021, GJI). Seismic\nattenuation represents the anelastic properties of rocks (Romanowicz & Mitchell, 2015, Treatise\non Geophysics; Takei, 2017, AREPS). Therefore, the seismic attenuation map determines the\nspatial pattern of effective viscosity (Fig. R1) but leaves the absolute value undetermined. In\npractice, we followed the same approach in Liu & Hasterock (2018, Science) to convert seismic\nattenuation to effective viscosity. The best-fitting conversion coefficients were determined by\nmatching the predicted maximum horizontal compression (S ) directions to the observation\nHmax\n(Levandowski et al., 2018, Nature Geoscience; Lund Snee & Zoback, 2020, Nature\nCommunications). In the revised manuscript, we have deleted all statements about the impact of\nthe effective viscosity structure on crustal stress to avoid circular arguments.\nSurface deformation patterns, such as the GPS-measured surface motion and deformation rate,\nare the goal of our prediction. Therefore, the observed surface deformation patterns were not\nused to construct the effective viscosity structure. Although most previous studies took GPS\nor geological measurements as known surface kinematic conditions to calculate depth-integ\nrated lithospheric effective viscosity structure (e.g., Flesch et al., 2000, Science; Ghosh an\nd Holt, 2012, Science; Gosh et al., 2013, JGR; Bahadori et al., 2021, Nature Communications),\nGhosh et al. (2013, JGR) showed that this approach could not fully reproduce the input\nGPS measurements. This observation calls for an alternative approach that is independent of G\nPS measurements to construct the three-dimensional (3D) effective viscosity structure, such a\ns our approach. We also tested the sensitivity of predicted surface motion to the effective visco\nsity structure, and it only shows weak sensitivity in the plausible regions. Therefore, the inclu\nsion or exclusion of geodetic measurements in the construction process does not impact the be\nst-fitting effective viscosity structure.\nIn practice, the long crustal weak zone was revealed and converted from the seismic attenuation\nmap showing high attenuation in this region (Fig. R1). We have revised line 520 (lines 1224-1225\nin the revised manuscript) to clarify that this crustal weak zone is in the input viscosity structure,\nwhich is converted from seismic attenuation.\nThe wording in the original “Methods” section and some places in the main text are confusing.\nTherefore, we have substantially revised the “Methods” section and the main text to clarify how\nwe construct the lithospheric effective viscosity structure and deleted all statements about the\nimpacts of the effective viscosity structure on crustal stress.\nIn conclusion, the takeaway of this paper is important, but the strength of their article is not clear\nin its current form due to the lack of quantitative comparisons between the different models (and\ntherefore, plate driving forces, which they discuss quite a bit in the Introduction to set up their\nproblem). Also, it is extremely difficult to keep track of all the models and their different\nassumptions, inputs and results.\nLastly, and perhaps most importantly, the lack of clarity in writing and terminology will prevent this\nmanuscript from having the impact that it likely could have given the importance of the result. My\nrecommendation is that a significant amount of manuscript revision needs to be done to make the paper\nmore intelligible, as it is extremely difficult to read and understand in its current form.\nResponses: We really appreciate your constructive insights. We have substantially revised the\nmanuscript for greater clarity and added more quantitative evaluations of our model prediction to\nbetter support our conclusions. Particularly, we have added quantitative comparisons between our\nmodel prediction and observation to show the effects of different parameters tested in this study,\nand the parameters and results are listed in the supplementary table S1. In addition, we have\nintroduced a quantitative evaluation of the relative importance of three different proposed driving\nforces in driving crustal deformation (Fig. 9). In the revised supplementary table S1, we have also\nincluded detailed setups of all models. We hope the revision has a better logical flow and can easily\nbe assessed.\nDetailed Comments:\nAbstract\nL21: “Particularly, the interaction between the mantle flow and lithospheric thickness step along\nthe eastern boundary of Basin and Range represents a key driving mechanism for localized\nintracontinental deformation.” Given that the motivation in the Intro is revealing the contributions\nof plate boundary forces, lithospheric body forces, and viscous drag to deformation, which\ncategory does this fall under? I could see this being a lithospheric body force (normal force exerted\non LAB topography from flow) and/or a viscous drag. It’s not clear to me that this study actually\nilluminates the contributions between the stated body forces. I think the motivation should be\nchanged to clearly apply to the conclusions, or the conclusions should be put in the context of the\nmotivation.\nResponses: Thank you for bringing this up! We have changed the abstract and introduction part to\nshow that our motivation is exploring the impacts of complex 3D lithospheric structures on crustal\ndeformation and stress state, as it is an unexplored but important aspect. The interaction between\nthe mantle flow and lithospheric thickness step will primarily change the contribution of mantle\nconvection to crustal deformation by changing the magnitude of crustal stress caused by the plate\nboundary effect and basal tractions.\nL26: uplifting -> uplift\nResponses: We have revised this sentence and the word “uplifting” has been deleted.\nL26-28: The first part of this sentence is grammatically unclear. I assume you mean that geological\nand geophysical evidence show clear signs of extension.\nResponses: Yes, here we mean the geological evidence shows clear signs of long-term extension.\nWe have revised this sentence for greater clarity (lines 48-50).\nL33-34: This sentence needs to be linked clearly to the stated 3 categories of driving forces.\nAlternatively, you could walk back the focus on the 3 categories, and simply state that you are\nexploring something that is largely unexplored in the context of these numerical models.\nResponses: Thanks for this suggestion. We have revised this sentence to show that the complex\n3D lithospheric structure affects all three proposed driving forces, and its effects were unexplored\nin previous studies (lines 55-57).\nL42-44: There has been a lot of research done on the controlling forces of this change in plate\nmotion and is largely attributed to plate boundary forces. An example paper exploring these forces\nshould be cited (perhaps a geodesy paper).\nResponse: We have revised this paragraph to acknowledge previous geodetic studies and a citation\nwas added in line 83.\nFig. 1: I recommend adding the boundaries of the Northern, Central, and Southern Basin and Range\nprovinces since you discuss them later in the paper. You should also cite the reference for LAB\ndepth.\nResponse: Thank you for the suggestions! We have added the boundaries of the northern, central,\nand southern Basin and Range (B&R) based on Jones et al. (1992, Tectonophysics). The LAB\ndepth is a hybrid model constructed upon LITHO 1.0 and a body-wave tomography model\n(Schmandt & Lin, 2014, GRL). We have added details about how we constructed this LAB depth\nmodel in the revised “Methods” section (lines 1737-1744).\nL68: Need to clarify that PNW volcanoes are not intraplate. The main edifices are typical calc-\nalkaline subduction-related volcanoes and shouldn't be included with intraplate volcanoes.\nResponse: This is a great point, and we really appreciate it. In the revised manuscript, we have\ndeleted all volcanism-related content according to another reviewer’s suggestion.\nL83-84 and 112: Could you clarify what you mean by \"one-way\" coupling? Seems like a\ncontradiction. Does \"one-way\" coupling mean all forces are drag from asth to lithosphere? What\neffect does \"two-way\" coupling actually have? It seems that all it implies is that it allows for more\nspatial complexity in your model.\nResponse: This is a good question! The one-way coupling approach only allows the convecting\nmantle to exert horizontal shear along the flat LAB. This horizontal shear is usually calculated in\na global mantle convection model only considering long-wavelength deep density anomalies.\nNotably, the lithosphere, which has a laterally varying thickness in reality, does not influence the\nmantel flow in this approach. The scenario modeled by one-way coupling may not be realistic in\nregions with sharp lateral gradients in lithospheric thickness, such as the WUS. Besides, our recent\nhigh-resolution mantle convection model shows that lateral variations in lithospheric thickness\nprominently affect the shallow mantle flow beneath the WUS, and such interaction is critical to\nreproduce the observed flow-induced SKS splitting pattern (Zhou & Liu, 2018, EPSL). However,\nthis lithospheric effect on mantle flow is neglected in the studies using the one-way coupling\napproach. We have revised this part to point out the meaning and defect of one-way coupling in\nlines 120-126.\nIn our fully dynamic model, the lithosphere and convecting mantle are allowed to interact naturally.\nFor example, the lithospheric thickness step along the eastern boundary of B&R blocks the shallow\nmantle flow in our model (Fig. 3). Besides, the shallow mantle flow pushes the thick cratonic\nlithosphere through normal stress along the lithospheric thickness step (Fig. 7). This natural\ncoupling allows more spatial complexity in our model, and it is more realistic.\nL84, L97"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Summary\nThis manuscript is modified from one that I reviewed previously.\nRecommendation\nThe authors have addressed most of the issues that I raised, which led to improvements to the\nmanuscript. I would consent to the publication of the paper provided that they adequately address\nthe issues detailed below.\nDetailed Review:\nDespite this version’s improved organization, I still find it difficult to understand the relationship\nbetween the many models that the authors present by only reading the main text. I thus suggest\nmoving Table S1 to the main text because it can greatly improve the clarity of the paper. Besides,\nin their response, they claim to have changed “Empty” in Table S1 to “None” following my\nsuggestions, but this is not true. “None” or “Not included” are more appropriate than “Empty”.\nHow the maps in Fig. 9 are computed are unclear to me. Lines 499–501 gives a brief description,\nbut the sentence is very difficult to understand. It says that they computed the ratios of the\nmagnitudes of the three types of driving forces, and the magnitudes are equivalent to the second\ninvariant of the stress tensor. My first question is: why the magnitudes are equivalent to the stress\ninvariants? Second, how did the authors separate the contributions of the three driving forces to\nthe stress field? Third, against what reference values are the ratios computed?\nI don’t see the point of presenting the model M10. This issue contributes to my general\nimpression that the paper lacks a clear focus: It seems to me merely a laundry list of all findings,\nbig and small, from an exhaustive parameter study, whereas I believe that papers suitable for\nhigh-impact journals such as Nature Communications should tell a cohesive story with a clear\nfocus. Both the other reviewer and I have alluded to this in our reviews, but the authors seem to\nbe only willing to make limited changes. I personally find this frustrating and strongly suggest\nthat the authors further trim the paper by removing redundant materials.\nLines 454–475: I also find this paragraph redundant because the focus of the paper is, at least\nfrom my reading it, how the complex lithospheric structures in the WUS control the intraplate\ndeformation in the region, NOT the origins of these complex structures. I thus suggest deleting\nthis paragraph.\nAlthough the authors claim that they have fixed all grammatical errors and typos, I still found a\nfew, which are listed below. There are likely more of them, so I urge the authors to carefully\nproofread their manuscript again.\nLine 213: “implying for” is grammatically incorrect. “for” should be deleted.\nLine 360: Should “Fig. 5e” be “Fig. 5h” instead?\nLine 363: It’s better to replace “considering” with “adding…to”.\nLine 388: “incorporates” is more appropriate than “considers”.\nLine 373: What “mantle convection” means here is unclear. The difference between M6 and M9\nand the other models seem to be the heterogeneous lithospheric structures, so “heterogeneous\nlithospheric structures” may be more appropriate here.\nLines 514–515: The “and” between “solutions” and “geophysically” should be replaced with a\ncomma.\nFig. 5: Why does M9 seem to have worse GPS misfits than M6? Why do Southern California\nand the area near the Mendocino Triple Junction have significantly higher GPS misfits\nconsistently in all scenarios?\nThe authors use the expression “X MPa” in the legends of a few figures, which looks very much\nlike a typo at first glance. I suggest adding a legend for each subplot instead.\nM3 should be labeled in Figs. 7 and 8 so that readers can easily see what models they are. Why is\n“with” italic in the caption of Fig. 7?\nFig. 8b: The bars are difficult to see. Consider increasing their thickness.", "reviewer_2": "Overall, I think the authors have done well in addressing my major concerns on the clarity of the\nmanuscript. I very much appreciate their statistical quantification of the results, and still strongly\nsuggest they find a way to implement this table in some form into their manuscript (perhaps just\nadding the quantitative measures of fit to some portion of the model result figures would be\nsufficient). I also suggest they add a Conclusion section where they concisely wrap up their major\nconclusions from their thorough search through the importance of the different tested parameters.\nI have some minor comments below, but find this manuscript suitable for publication once these\nare addressed.\nBelow are responses to the researchers' corrections, with original comments in black, their\nresponse in blue, and my follow-up comments in red. Comments that are deleted from prior review\nhave been sufficiently addressed\nL33-34: This sentence needs to be linked clearly to the stated 3 categories of driving forces.\nAlternatively, you could walk back the focus on the 3 categories, and simply state that you are\nexploring something that is largely unexplored in the context of these numerical models.\nResponses: Thanks for this suggestion. We have revised this sentence to show that the complex\n3D lithospheric structure affects all three proposed driving forces, and its effects were unexplored\nin previous studies (lines 55-57).\nL34: Should change \"affects\" to \"is linked to\"\nL42-44: There has been a lot of research done on the controlling forces of this change in plate\nmotion and is largely attributed to plate boundary forces. An example paper exploring these forces\nshould be cited (perhaps a geodesy paper).", "reviewer_author_response": "Response letter\nWe greatly appreciate the time and effort the reviewers devoted to providing valuable and\ninsightful feedback on our work. The constructive comments and suggestions are beneficial in\nenhancing the quality of our manuscript. In response to the feedback, we have tried our best to\nincorporate all the suggestions and comments, with tracked changes in the revised manuscript.\nOnce again, we sincerely appreciate the reviewers’ constructive inputs to our work. We are\ncommitted to addressing any further concerns the reviewers may have.\nBelow are our detailed point-to-point responses to the reviewers’ comments and suggestions, with\noriginal comments in black and our response in blue.\nTo reviewer #1:\nThe authors have addressed most of the issues that I raised, which led to improvements to the\nmanuscript. I would consent to the publication of the paper provided that they adequately address\nthe issues detailed below."}, "author_response": ""}, {"round": 3, "reviewer_comments": {"reviewer_1": "Recommendation\nI was not expecting to review this manuscript again because in the last round of review, I said I\nwould consent to the publication of the manuscript provided that the authors address my\ncomments. Nonetheless, the editor insisted on asking me to review it for a third time, suggesting\nthat she also has concerns over the work. I confess that I was very reluctant to do so because I\nhad already spent so much time on it, and that the experience of reading the manuscript is far\nfrom being pleasant, mostly due to its poor organization, about which I have complaint many\ntimes. Nonetheless, out of the sense of responsibility to the community, I decided to take up the\nunpleasant task once again. However, I was soon shocked by how disrespectful the authors are to\nmy time: Table 1, which was brought to the main text upon my request and is instrumental in\nclarifying their bewildering parameter-test procedure, contains two IDENTICAL rows (5 and\n10)! I read the two rows multiple times to make sure that I didn’t hallucinate because I simply\ncould not believe that such an obvious typo could appear in such an important place in a\nmanuscript that has already gone through two rounds of review. This incident, together with the\nmany typos and grammatical errors that I caught in the previous two rounds, suggests to me that\nthe authors are extremely negligent, which inevitably leads one to doubt the reliability of their\nwork. I am also furious about their complete lack of respect for my precious time spent in\npainstakingly reviewing their work. I thus refuse to review the manuscript again unless I receive\nthe guarantee that my time spent in doing so will be respected."}, "author_response": "Response:\nWe sincerely apologize for the remaining problem in the manuscript and thank you for\npointing this out. In the last revision, we reduced the number of models presented in the main text\nand changed the order of models, which led to the remaining problem. In this revision, we have\ncarefully examined all text and display items and corrected the problem in Table 1.\nWe truly appreciate your precious time and effort in providing constructive comments and\nsuggestions to improve our manuscript’s quality. Again, we sincerely apologize for all the typos\nand grammatical issues in previous versions of our manuscript, and we express our most sincere\ngratitude for your time and effort devoted to helping us improve the manuscript."}]

Decision letter