Total Variation Regularization of Geodetically and Geologically Constrained Block Models for the Western United States

Geodetic observations of interseismic deformation in the Western United States provide constraints on microplate rotations, earthquake cycle processes, and slip partitioning across the Pacific–North America Plate boundary. These measurements may be interpreted using block models, in which the upper crust is divided into microplates bounded by faults that accumulate strain in a first-order approximation of earthquake cycle processes. The number and geometry of microplates are typically defined with boundaries representing a limited subset of the large number of potentially seismogenic faults. An alternative approach is to include a large number of potentially active faults bounding a dense array of microplates, and then algorithmically estimate the boundaries at which strain is localized. This approach is possible through the application of a total variation regularization (TVR) optimization algorithm, which simultaneously minimizes the L 2 norm of data residuals and the L 1 norm of the variation in the differential block motions. Applied to 3-D spherical block models, the TVR algorithm can be used to reduce the total variation between estimated rotation vectors, effectively grouping microplates that rotate together as larger blocks, and localizing fault slip on the boundaries of these larger block clusters. Here we develop a block model comprised of 137 microplates derived from published fault maps, and apply the TVR algorithm to identify the kinemati-cally most important faults in the western United States. This approach reveals that of the 137 microplates considered, only 30 unique blocks are required to approximate deformation in the western United States at a residual level of <2 mm yr −1 .


INTRODUCTION
Understanding the geometric complexity of active fault systems has been an extant question since the advent of plate tectonics and is necessary for both quantitative seismic hazard assessment and understanding continental tectonics. Quantitative models of plate kinematics at global scales have included 12 (DeMets et al. 1990) to 19 (Sella et al. 2002) kinematically distinct plates, and phenomenological plate identifications have suggested at least 52 (Bird 2003). The primary challenge in developing global models that incorporate both plate-boundary scale (>1000 km) and regional scale (<1000 km) microplate and fault system geometries is uncertainty regarding the level of discretization necessary to explain the kinematics of active fault zones. Over the last 30 yr, Global Positioning System (GPS) observations have provided a quantitative description of nominally interseismic continental deformation with which to address this challenge and constrain plate boundary kinematics. However, elastic strain accumulation across locked faults produces smooth GPS velocity gradients (e.g. Savage & Burford 1973), * Now at: U.S. Geological Survey, Menlo Park, CA 94025, USA. making it potentially challenging to identify the slip contribution of individual faults within a fault system without explicitly defining fault geometries and modelling interseismic earthquake cycle effects.
To address this potential ambiguity in model realizations of fault system geometry, kinematic descriptions of western United States tectonics were developed satisfying path integral constraints and slip rate observations from geology (Weldon & Humphreys 1986;Minster & Jordan 1987) and geodesy (Minster & Jordan 1987) to quantify the potential influence of poorly documented deformation sources. This kinematic consistency constraint has formed the foundation of block models (Matsu'ura et al. 1986;Bennett et al. 1996Bennett et al. , 1997Souter 1998;McClusky et al. 2001;Murray & Segall 2001;Becker et al. 2005;McCaffrey 2005;McCaffrey et al. 2007;Meade & Hager 2005;Johnson & Fukuda 2010;Spinler et al. 2010;Chuang & Johnson 2011;Hammond et al. 2011;, in which the crust is divided into a discrete number of microplates bounded by mapped faults and combined with a first order approximation of quasi-static earthquake cycle deformation. Assessing fault system complexity in a block model has often been approached through inspection (e.g. Meade & Hager 2005), 714 E.L. Evans, J.P. Loveless and B.J. Meade limiting analysis to faults with known geological fault slip rate estimates (Chuang & Johnson 2011), or by considering a small number of geometrically variable block model configurations (e.g. d'Alessio et al. 2005;Spinler et al. 2010). More recently, clustering algorithms have been applied to interseismic GPS velocities in the San Francisco Bay Area (Simpson & Thatcher 2012) and the eastern Mojave (Savage & Simpson 2013) to identify statistically distinct regions of block-like motion. Cluster analysis provides an algorithmic approach to block identification though it neglects earthquake cycle processes and is designed for cases where GPS observation coordinates are distant from the location of the Euler pole of the block on which they lie.
Here we develop block models containing a dense array of microplates with boundaries based on mapped fault traces in the western United States and apply a total variation regularization (TVR) optimization algorithm as a tool for block model selection. This quantized state vector estimation technique enables grouping of adjacent small microplates into fewer, larger blocks and identifies the kinematically most important block boundaries. In this way, we seek to algorithmically determine the simplest plate boundary block model (i.e. with the fewest blocks) that fits the data at a given level of resolution.

BLOCK MODELLING AND Q UA N T I Z AT I O N W I T H T O TA L VA R I AT I O N R E G U L A R I Z AT I O N
In the block model formulation of Meade & Loveless (2009), interseismic velocities are represented as a linear combination of block rotations and elastic strain accumulation due to slip deficit on block-bounding faults and spatially variable strain accumulation on triangular dislocation elements (TDEs, Meade 2007), where v I are the interseismic velocities, v B are the velocities due to block rotation, v E are velocities due to a first order approximation of elastic earthquake cycle effects on locked faults and due to elastic strain accumulation on TDEs and v R are residual velocities due to model error and observational noise. Written in terms of block rotations and elastic strain accumulation, the forward problem is: where G B contains the partial derivatives functions of velocity observations with respect to the block rotation vectors, G SD contains the integrated effects of Green's functions associated with slip deficits on block-bounding faults, G t contains Green's functions for slip deficit on triangular dislocation elements, is the vector of block rotations and t contains the slip deficit values on partially coupled triangular dislocation elements. In this formulation, all block-bounding faults, except those explicitly represented with triangular dislocation elements, are assumed to be fully locked and accumulating strain at the full estimated fault slip rate. For the commonly considered block model problem, the block motions and TDE strain accumulation rates are considered unknown and solved for in a least squares sense, parameter β, m est = [ est , t est ] ′ and d = v I are interseismic velocities from annual to decadal scale GPS observations (Meade & Loveless 2009). The block models used here are constrained by horizontal motions only because a kinematically consistent block model formulation is not currently extended to vertical deformation.
Within a system of rotating microplates, a fault bounding two blocks will be active (i.e. have a non-zero slip rate) when the blocks have different rotation vectors, whereas if two blocks have an identical Euler pole location and rotation rate, those two blocks effectively behave as a single larger block, and the slip rate on the bounding fault will go to zero (Fig. 1). We may wish to find a solution in which many blocks have identical rotation vectors, that is, in which is quantized. Estimating quantized rotation vector solutions provides a means for approaching ambiguity in block model selection: we may include all possible faults in a dense array of microplates, and then algorithmically estimate the boundaries at which strain is localized based on geodetic observations of interseismic deformation.
It is possible to estimate quantized solutions directly with a regularization technique developed for edge sharpening in image processing: TVR (Rudin et al. 1992;Chambolle 2004). TVR simultaneously minimizes the L 2 norm of the data residuals and L 1 norm of the variation in the estimated state vector. Applied to 3-D spherical block models, TVR minimizes the L 1 norm of variation in the block rotation vectors, therefore localizing fault slip on the boundaries of these larger blocks, where λ controls the strength of the regularization term, and D is a linear differential operator. Construction of matrix D is discussed at the end of this section. Absolute value, L 1 , regularization methods designed to recover sparse state vectors have been applied in reflection seismology for decades (e.g. Claerbout & Muir 1973;Santosa & Symes 1986;Yao et al. 2011), and recently applied to geodetic imaging of sharp coseismic slip distributions . Under many circumstances, the L 1 norm likely recovers the L 0 pseudo-norm, which gives the number of non-zero elements in the vector of interest (Candès et al. 2006;Donoho 2006). Applied to earthquake slip distributions, L 1 regularization minimizes the number of non-zero elements in the solution vector, and produces a compact representation of slip . Applied to block models, applying the L 1 norm to the vector Dm limits the number of nonzero differences in the solution vector m, resulting in a grouped or quantized solution. The L 1 regularization is non-linear and a global minimum can be found using convex optimization methods (Boyd &  Vandenberghe 2004). Several solvers exist for variations of L 1 regularization including TVR. We have modified a one-dimensional version of the TVreg package for Matlab (Boyd & Vandenberghe 2004;Jensen et al. 2012; http://www.imm.dtu.dk/∼pcha/Regutools/TVreg.m) to accommodate arbitrary problem sizes and constraint matrix structure.
In terms of a block model, every three rows of D correspond to a block constraint, and each column corresponds to a rotation vector component of a block: D has dimensions of three times the number of constraints by three times the number of blocks. Each constraint limits differences in the rotation vector between two adjacent blocks. For example, if block i is adjacent to block j, we penalize models where the difference in the rotation vectors between these blocks, i − j is not zero. For each block, we calculate the difference between all geographically adjacent blocks by specifying separate rows in the difference matrix for each of the three components of the rotation vectors. For example, the constraint matrix, C, for single pair of neighbouring blocks can be written as, Each of the constraint matrices, C, for adjacent block pairs, are stacked to form difference matrix, D.

GPS DATA, GEOLOGICAL SLIP R AT E S , A N D B L O C K M O D E L C O N S T RU C T I O N
We constrain deformation in the western United States using horizontal interseismic velocities at 1686 locations (McClusky et al. 2001;Shen et al. 2003;Hammond & Thatcher 2005;Williams et al. 2006;McCaffrey et al. 2007;; and the Plate Boundary Observatory Network velocity field, http://pboweb.unavco.org; Fig. 2a). These velocity fields are combined into a common reference frame by minimizing velocity misfit between collocated stations in the fields, using a six-parameter (three-component rotation and three-component translation) transformation , GSA Data Repository item 2011305). As discussed in Section 2, we focus our analysis on horizontal GPS velocities. We have omitted velocities north of the Canadian border because we do not model potentially complex deformation in western Canada or along the Aleutian trench, and velocities near Yellowstone that contain contributions from nontectonic deformation.
We additionally constrain western United States tectonics with 48 geological slip rate estimates (Fig. 2b, Table S1). These observations are selected from the catalogue compiled by Dawson & Weldon (2013), and we include three slip rates on faults outside of California on the Annie Spring fault (Bacon et al. 1999), the Petrified Springs fault (Wesnousky 2005) and the Wasatch fault (Jewell & Bruhn 2013). We have considered only geological rates >1 mm yr −1 . To additionally account for epistemic uncertainties in the geological rates that may not be incorporated in the reported rates (e.g. Bird 2007;Gold et al. 2009;Zechar & Frankel 2009), we expand the uncertainty magnitudes to be at least 50 per cent each respective reported rate. Where a single fault has multiple geological slip rate estimates associated with it, we select a single published rate to avoid ambiguities in combining palaeoseismic uncertainties (e.g. McQuarrie & Wernicke 2005). The rates, reported uncertainties, applied uncertainties, and sources for the geological slip rates are provided in Table S1.
In the block model, we incorporate geological slip rates as a priori slip rate observations (Meade & Loveless 2009), which modifies the components of eq. (3) such that: 716 E.L. Evans, J.P. Loveless and B.J. Meade  where G ap contains the partial derivative functions of the a priori rates with respect to the block rotation vectors, and vector s ap contains the a priori geological slip rate observations. In this modification, the data-weighting matrix W contains additional diagonal terms describing the relative weighing of the geological slip rates relative to the GPS velocities (Meade & Loveless 2009, eq. 13). We construct W such that the total influence of the geological slip observations are down-weighted relative to the geodetic velocity observations by a factor of 10.
We develop a block model comprised of 137 microplates based on Quaternary fault maps (Jennings 1994;U.S. Geological Survey et al. 2006;Figs 2b and 3). Similar to other block modelling studies, all locking depths on block bounding faults are assumed to be 15 km (e.g. Meade & Hager 2005;Johnson 2013). We have incorporated two exceptions in which we model creeping faults with a locking depth of 0 km: on the Hayward fault in the east San Francisco Bay Area and on the creeping segment of the central San Andreas fault.
We divide the block modelling procedure into two steps: (1) we use TVR to reduce the number of uniquely rotating blocks based on the GPS observations and geological slip rate constraints. The grouped blocks define the block geometry to be used in modelling deformation in the western United States. Faults bounding these larger blocks are considered active, and the zero-slip-rate interior faults are removed from the block model geometry. In this step, we use a modified velocity field with a previously estimated elastic contribution due to Cascadia subduction ) removed, so that the resulting velocity field may be attributed exclusively to block rotations and earthquake cycle processes on block-bounding faults.
In step (2), we use the block geometry selected by step (1) for a full block model solution constrained by the full velocity field and geological slip rate constrains on the selected fault geometry with . Flowchart outlining workflow for TVDN model selection and block model estimation. We use TVDN as a method of block geometry selection based on original dense block geometry and a modified velocity field (modified to remove the elastic signal due to coupling on the Cascadia subduction zone). The selected block geometry is then fed into the full block model using the full velocity field to estimate block rotation, fault slip rates, and coupling on the subduction zone. a traditional block model (eq. 1) using a weighted least squares estimator (eq. 3) to estimate block rotations, fault slip rates, and spatially variable slip on triangular dislocation elements representing the Cascadia subduction zone (McCrory et al. 2009;Fig. 2b). For simplicity, we identify a single smoothing parameter β = 1000 for the subduction zone, which we apply to all models considered here, although this may be varied as well to obtain desired smoothness and goodness of fit properties of the subduction zone solution.
For the full block model inversion we assume a priori rotation constraints on the motion of the Juan de Fuca plate due to a lack of direct GPS observations [Euler pole location −26.6 ± 0.1 • N, 69.6 ± 0.1 • E, and rotation rate 0.810 ± 0.001 deg Myr −1 (Miller et al. 2001)]. This procedure is outlined in the flowchart in Fig. 4.

QUANTIZED BLOCK MODEL R E S U LT S
Total variation regularization reduces the effective number of blocks directly by simultaneously estimating a set of block rotations that are quantized. The parameter λ controls the strength of the regularization such that increasing λ reduces the number of unique rotation vectors (i.e. the number of kinematically distinct blocks). Very small values of λ do not reduce the total number of independently rotating blocks. For large values of λ, model fit to the GPS observations is outweighed by the total variation constraint, resulting in models with very few independently rotating blocks, and therefore very few active faults. We consider block model solutions for λ ranging from λ = 50 to 3150, which incrementally reduces the total number of blocks in the plate boundary from 137 to 3, while data misfit, quantified in terms of mean residual velocity (MRV), increases from 1.4 to 4.1 mm yr −1 (Fig. 5). We also track regional MRV in southern California, the San Francisco Bay Area, the Basin and Range and the Pacific Northwest. Vertical dashed lines in Fig. 5 identify λ values of models that we discuss in this section.
As a benchmark and end-member, we first consider the λ = 3150 model, in which the plate boundary is made up of only three blocks: a Pacific block, a North America block, and a Juan de Fuca block (Fig. 6). In this extreme case, the best-fitting boundary between the Pacific and North American blocks follows the San Andreas fault for much of its mapped trace between the Cerro Prieto fault in Mexico and the central San Andreas fault, with one exception. South of 34.2 • N in southern California, slip follows the San Jacinto fault west of the San Andreas fault, reconnecting with the Imperial fault south of 33 • N. North of the central San Andreas As λ increases, the number of unique Euler poles (N) decreases (blue line), and the mean residual velocity (MRV) increases (dark red line). Regional MRV for southern California (light red dashed line) follows the overall plate boundary MRV, while the regional MRV for the San Francisco Bay Area (light red solid line) is consistently larger. Dashed black lines identify λ values of models that we discuss in Section 4: λ = 650, reference model λ = 850 and λ = 1150. Regional mean residual velocities are calculated on the area shown in Figs 9-15. fault, slip on the boundary steps east onto the Valley Margin fault, and follows the Greenville and Bartlett Springs faults, eventually merging with the Cascadia trench offshore of southern Oregon ∼400 km north of the Mendocino triple junction. Slip rates on this end-member plate boundary are 28-33 mm yr −1 , and fit the data poorly, with a MRV = 4.1 mm yr −1 (Fig. S1). Regional misfits in Southern California, the San Francisco Bay Area, and the Basin and Range are higher, with local MRVs of 4.9, 6.5 and 4.2 mm yr −1 , respectively.
In this end-member case, slip on the Cascadia subduction zone accommodates strain that would otherwise be accommodated by faults in the northwest United States and in the Basin and Range. This is most apparent in a very high slip patch slipping at >40 mm yr −1 of normal-sense slip at the southern end of the Cascadia subduction zone (Fig. 6).

Reference and comparison models
We select a reference model of λ = 850 that contains 30 blocks, and explains GPS observations with a MRV of 1.7 mm yr −1 . Of the models considered, λ = 850 is the largest λ value for all regional MRVs are less than 2 mm yr −1 , with San Francisco Bay Area MRV = 1.8 mm yr −1 , Southern California MRV = 1.9 mm yr −1 , Basin and Range MRV = 1.2 mm yr −1 , Pacific Northwest MRV = 1.6 mm yr −1 and overall MRV = 1.7 mm yr −1 (Fig. 7). To demonstrate the range of possible fault geometries and slip rate estimates, we additionally consider a less quantized model with λ = 650, containing 41 blocks for an MRV of 1.6 mm yr −1 (Fig. 7a); and a more quantized model λ = 1150 (the largest λ value to fit GPS velocities with MRV < 2 mm yr −1 ) containing only 16 blocks, and with an MRV of 1.9 mm yr −1 (Fig. 7c).
In general, the λ = 650 model produces lower slip rates (on more faults) than the λ = 850 reference model, and λ = 1150 model produces higher slip rates (on fewer faults). The residual velocity fields for the reference and comparison models do not contain large-scale systematic orientations or magnitude (Fig. 8), with consistently smaller residuals in the λ = 650 model and the highest residuals in the λ = 1150 model. Among all the models, the highest residuals occur along the southern and central San Andreas fault (especially in the λ = 850 and 1150 models), northern California near the Oregon border, and Walker Lane (especially in the λ = 1150 model). Residuals are consistently low (<2 mm yr −1 ) adjacent to the Cascadia subduction zone.
The block model formulation produces uncertainty estimates on fault slip rates estimated formally for a given block geometry. These uncertainties therefore do not reflect uncertainties in block and segment geometry, and are therefore much smaller than differences in slip rate on a given fault between the models considered here. Because of this, when reporting model slip rates for a given block geometry, we do not report the formal uncertainties as reported by the block model.
For purposes of this discussion, we focus primarily on the strikeslip component of fault slip rates, especially on the westernmost edge of the plate boundary in California; and discuss specific examples of dip-slip and tensile deformation along the Cascadia subduction zone and in the Basin and Range.

Southern California
Between the reference and comparison models, we identify 7-16 rotating blocks in southern California, with 11 blocks in the λ = 850 reference model (Fig. 9). The reference model (λ = 850) fits GPS observations in southern California with a MRV of 1.7 mm yr −1 , compared with 1.7 mm yr −1 (λ = 650) and 2.1 mm yr −1 (λ = 1150). This survey of the block model selection results in southern California begins in the eastern California shear zone and progresses west to the offshore Borderlands faults (Fig. 10).
The     Harper fault zone, the Blackwater/Calico fault, and the Goldstone Lake fault. In the reference model, 12 mm yr −1 of right lateral slip is fed into the ECSZ by the Hidden Springs fault (Jennings 1994 (Fig. 9a).
The resulting block geometries in the reference and comparison models contain thin oblong north-south-oriented blocks that span the length of the ECSZ. Estimated strike-slip rates reach up to 13 mm yr −1 , consistent with previous geodetic studies (Gan et al. 2000;McClusky et al. 2001;Peltzer et al. 2001;Savage et al. 2001;Meade & Hager 2005), but higher than geological rates (e.g. Oskin & Iriondo 2004;Oskin et al. 2007;Oskin et al. 2008). At λ = 850, deformation in the ECSZ is localized on four of seven possible structures, and slip is not distributed evenly between them (Fig. 9b). Our results are consistent with previous results in which ECSZ slip is distributed across a small number of faults (McClusky et al. 2001;Meade & Hager 2005). In the λ = 650 reference model, only one additional fault is active to the east, suggesting that adding additional faults to approximate continuous deformation may not be necessary to describe ECSZ motion (e.g. Peltzer et al. 2001;Chuang & Johnson 2011;Fig. 9a). However unlike Peltzer et al. (2001), we do estimate more than a single fault in the eastern Mojave.
In the reference model, we estimate 1-5 mm yr −1 of right lateral slip on most of the Blackwater fault, increasing to 10 mm yr −1 in the south (Fig. 9b). This range is consistent with the ∼7 mm yr −1 of right lateral slip estimated by Peltzer et al. (2001). East of the Blackwater fault, we estimate 3-8 mm yr −1 of right lateral slip Goldstone fault, which is the easternmost active ECSZ fault in the reference model. Immediately west of the Blackwater fault, we estimate 5 mm yr −1 of right-lateral slip on the Harper fault zone. The next active faults to the east are the Lenwood/Lockhart faults, slipping right laterally at ∼6 mm yr −1 . We do not estimate any fault activity between the Lockhart and Lenwood faults and the southern San Andreas fault. Among the reference and comparison models, slip rates in the ECSZ vary within an order of magnitude, often with short-wavelength slip rate variations within a single model. This is apparent in all three models on Lenwood and Harper faults, which have high slip rates (>15 mm yr −1 ) in the north near the intersection with the Garlock fault but low slip rates (<10 mm yr −1 ) in the south (Fig. 9).
None of the models considered contain an active Garlock fault along its entire ∼200 km length (Fig. 9). Estimated slip rates on the Garlock fault are largely heterogeneous, reaching a local maximum near 118 • W of 12-18 mm yr −1 left lateral slip (17 mm yr −1 in the reference model). East of the intersection with the Blackwater fault, slip rates on the active segments of the Garlock fault do not exceed 3 mm yr −1 . Instead the reference and comparison models all contain a fully active left-lateral White Wolf fault (e.g. Bawden et al. 1997).
We estimate an active southern San Andreas fault along its entire trace with the fastest San Andreas rates on the Cerro Prieto fault in northern Mexico of 41 mm yr −1 . In the reference and λ = 650 models, 5-8 mm yr −1 of slip is transferred westward along the Elsinore fault. We estimate 9 mm yr −1 on the San Jacinto in the reference model, and up to 11 mm yr −1 in the λ = 1150 model (Fig. 9). In all three models slip follows Superstition Mountain section but not the Superstition Hills section. These San Jacinto rates are consistent with geodetic studies estimating 12 mm yr −1 (Lisowski et al. 1991;Lindsey et al. 2013), and considerably lower than the 25 mm yr −1 slip rate estimated by Lundgren et al. (2009) (who estimated ∼17 mm yr −1 on the San Andreas), and to Fialko's (2006) estimate of 19-20 mm yr −1 .
North of the juncture with the San Jacinto, slip on the Imperial segment of the San Andreas fault is 22-24 mm yr −1 , decreasing to 10-14 mm yr −1 along the Indio segment. The San Andreas fault exhibits a local minimum along the San Bernardino segment, with rates as low as 5 mm yr −1 in all three models. North of merging again with San Jacinto at 34.2 N, in the reference model we estimate 17 mm yr −1 along the Mojave segment and 31 mm yr −1 along the Carrizo segment (Fig. 9).
West of the San Andreas fault and offshore, the Agua Blanca fault in northern Baja California feeds as much as 12 mm yr −1 of slip into borderlands and offshore faults (Fig. 9). In the reference model, the Palos Verdes fault, and the Newport/Inglewood fault zone are active at 0.5-5 mm yr −1 . We estimate 4.2 mm yr −1 of right lateral slip on the Palos Verdes fault, which is consistent the geologically constrained rate of 4 mm yr −1 (Brankman & Shaw 2009). In the Los Angeles basin, the reference model contains an active Newport/Inglewood fault and Compton thrust with up to 6.5 mm yr −1 of convergence. This slip rate may be consistent the combined influence of actively uplifting Compton thrust (Leon et al. 2009), the Puente Hills blind thrust (Leon et al. 2007) and the Elysian Park thrust (Davis et al. 1989). In the λ = 650 model, slip is distributed across additional offshore faults: the San Diego trough fault zone and the San Clemente fault with right lateral slip rates of 2 and 4 mm yr −1 , respectively. In contrast, the λ = 1150 model has only one active structure in southern California west of the San Andreas fault consisting of the Rose Canyon-Palos Verdes faults with a strike-slip rate 4.4 mm yr −1 (Fig. 9c).
North of the Los Angeles basin, the reference model identifies 5 mm yr −1 of right lateral slip on the Malibu Coast fault, and 3 mm yr −1 of extension on the Ventura fault (Fig. 9b). Extension in the Ventura region is surprising given recent observations of uplift (Hubbard et al. 2014). We estimate localized shortening of up to 12 mm yr −1 on the San Gabriel fault, which extends north from the Sierra Madre fault. This shortening is more complicated in the λ = 650 model, in which both the southern San Gabriel fault 720 E.L. Evans, J.P. Loveless and B.J. Meade and the Sierra Madre fault are active with up to 20 mm yr −1 of oblique convergence on the San Gabriel fault immediately adjacent to the same magnitude of oblique extension on the Sierra Madre fault (Fig. 9a). These high rates and short wavelength variations are physically unlikely, suggesting that in this λ = 650 case, the total variation constraint is not strong enough to avoid over fitting the GPS data.
In the reference model at latitude 35-36 N • adjacent to the central San Andreas fault, the reference model and λ = 1150 models contain no active subparallel structures west of the San Andreas fault (Fig. 9). The λ = 650 model does contain active San Juan and Rinconada fault zones immediately to the east of the San Andreas, as well as an active offshore San Lucia Bank fault zone (slipping 8, 6 and 1 mm yr −1 , respectively; Fig. 9a). This fault geometry leads to high block rotations east of Parkfield, with extension rates on the west branch of the San Juan fault zone of up to 23 mm yr −1 immediately adjacent to convergence rates of 22 mm yr −1 on the east branch, again indicating over fitting in the λ = 650 model. North of the Carrizo segment, we estimate 31 mm yr −1 at Parkfield, consistent with previous geodetic estimates (e.g. Murray & Langbein 2006), and slightly faster than the geological estimate of 24.8 mm yr −1 (Toké et al. 2011). Slip splits north of Parkfield with 33 mm yr −1 on San Andreas fault and 3 mm yr −1 of left-lateral slip on a subparallel structure to the east, the Valley Margin fault .
Whereas the λ = 650 constraint appears to over fit GPS observations, the total variation constraint at λ = 1150 is likely too large to effectively reproduce the geodetic data in some areas. For example, none of the offshore faults are active, leading to high residual velocities offshore of >3 mm yr −1 (Fig. 8). While individual geometries and slip rates in the comparison models may not capture all southern California deformation, southern California tectonics appear to be best described by the rotation of more than 7 but fewer than 20 rotating blocks.

San Francisco Bay Area
In the reference model (λ = 850) and comparison models (λ = 650, λ = 1150), slip in the San Francisco Bay Area is partitioned across four to seven closely spaced active faults, with six major faults bounding five rotating blocks in the reference model (Fig. 11). The reference model fits the regional Bay Area GPS observations with MRV = 1.8 mm yr −1 . This is relative to a regional MRV = 1.8 (λ = 650) and MRV = 2.5 (λ = 1150) in the comparison models. Active faults include the San Andreas fault, the San Gregorio fault, the Hayward fault, the northern Calaveras and West Napa faults, and the Greenville fault. East of the southern Greenville fault, we include an active structure that was also included in the model of d' Alessio et al. (2005) as the Valley Margin fault (Fig. 12). Active structures are consistent between the reference and comparison models, although the λ = 650 model contains an additional fault east of the Greenville fault, and the λ = 1150 model does not include an active peninsula segment of the San Andreas fault (Fig. 11).
There is considerable variation in slip rates between the reference and comparison models, in particular in the western Bay Area. In the comparison models, the Zayante fault transfers 3.5-10 mm yr −1 of right lateral slip from the central San Andreas fault across the Santa Cruz mountains to the San Gregorio fault in the southwest Bay Area (Figs 11a and c). The reference model, however, estimates left-lateral slip rates on the Zayante fault of 5 mm yr −1 , which then connects to an 11 mm yr −1 left-lateral San Gregorio fault (Fig. 11b). Adjacent to the San Gregorio fault in the reference model, we estimate 23 mm yr −1 of right lateral slip on the peninsula segment of the San Andreas fault, but this segment is not active at all in the λ = 1150 model. Previous block models of the San Francisco Bay Area have not included the Zayante fault , and its presence in the comparison models may be interpreted as consistent with active interseismic strain accumulation near the epicentre of the 1989 Loma Prieta earthquake (e.g. Shaw et al. 1994), although in this study its behaviour appears to be poorly constrained.
In the East Bay, we allow the Hayward fault to creep at its full estimated rate in all of the models by setting its locking depth to zero. Our slip rates estimated on the Hayward fault are 6-9 mm yr −1 , with 6.2 mm yr −1 in the reference model, similar to estimates from studies allowing partial creep (e.g. Bürgmann et al. 2000;d'Alessio et al. 2005;Schmidt et al. 2005;. In the reference model, the southern Calaveras fault transfers 9-10 mm yr −1 of slip from the creeping segment of the San Andreas to the southern Hayward fault, which then splits between the 6 mm yr −1 Hayward fault and a 5 mm yr −1 fault cutting across the south bay to reconnect with the peninsula segment of the San Andreas fault (Fig. 11b). The Hayward fault extends to the Rodgers Creek fault to the north, slipping up to 17 mm yr −1 , and to the Bloomfield fault slipping 4 mm yr −1 . The Bloomfield fault connects the Hayward to the northern San Andreas fault. However, fault behaviour in the northern San Francisco Bay Area is variable among the models considered. The Bloomfield fault is inactive in the λ = 650 model, but active in the reference and λ = 1150 models, in which the Hayward-Bloomfield-Rodger's Creek triple-junction accommodates up to 30 mm yr −1 of extension and compression on adjacent fault segments (Fig. 11).
East of the Hayward fault, the southern Greenville fault transfers slip from the central San Andreas fault, linking up with the northern Calaveras. In the reference model and λ = 650 model, the Valley Margin fault slips left-laterally at ∼5 mm yr −1 eventually connect-ing with the Greenville fault, although this fault remains right lateral in the λ = 1150 model (Fig. 11). At the intersection between the Greenville and northern Calaveras faults, rather than transferring slip directly, the Greenville and northern Calaveras faults bound a small a small block in the easternmost Bay Area (Fig. 11) in all three of the models considered. We estimate a slip rate on the northern Calaveras fault of 10 mm yr −1 , faster than the geological slip rate estimate of 5.6 mm yr −1 (Simpson et al. 1999), and geodetic estimates of 6.2-9.0 mm yr −1 (d' Alessio et al. 2005;. Calaveras slip rate is transferred to the West Napa fault to its north, which slips 8-10 mm yr −1 in all of the models (10 mm yr −1 in the λ = 850 reference model), and was the location of the recent M W = 6.0 South Napa earthquake (U.S. Geological Survey, earthquake.usgs.gov).
While there is considerable variability in fault geometry and slip rates, the reference and comparison models appear to constrain complexity of the Bay Area fault system as best described by 4-7 rotating blocks.

Basin and Range and Walker Lane
Slip is partitioned through the Basin and Range province similarly among the reference and comparison models. In the reference model (λ = 850), slip is bounded on the west side of the province by several faults in the eastern California shear zone and up through the Walker Lane belt (Stewart 1988), and on the east by the Wasatch fault (Figs 13 and 14). A fault following the Central Nevada seismic belt (Wallace 1984a) remains active in the western portion of the interior of the province, continuing up into Idaho. The λ = 650 model estimates similar slip partitioning, with some geometry differences in Walker Lane, and the λ = 1150 model does not include activity along the Central Nevada seismic belt (Fig. 13).
This set of active structures is similar to those described by Bennett et al. (2003), although we do not identify active slip on the Intermountain Seismic Belt (Smith & Sbar 1974). In the reference and comparison models, no slip rates on structures in the interior of the Basin and Range exceed 3 mm yr −1 , similar to slow rates estimated by Niemi et al. (2004;Fig. 13). In the reference model, ∼10 mm yr −1 of right-lateral slip is transferred from the eastern California shear zone through the eastern Walker Lane fault zone, with slip rates up to 8 mm yr −1 in Walker Lane, somewhat higher than of total shear of ∼7 mm yr −1 (Hammond & Thatcher 2007). The more quantized comparison model of λ = 1150 results in little deformation within Walker Lane, with 7 mm yr −1 of slip transferred instead from the eastern California shear zone into the Central Valley of California. Another 1-2 mm yr −1 of slip is accommodated on a structure in eastern Walker Lane (Fig. 13c).
At the eastern edge of the Basin and Range, the Wasatch fault is active in the reference and comparison models, and accommodates only 2-4 mm yr −1 of extension (2.8 mm yr −1 in the reference model). The interior structure along the Central Nevada Seismic belt accommodates a similar magnitude of extension: 0.5 mm yr −1 near the Idaho border to 3 mm yr −1 in western Nevada and 2 mm yr −1 of strike slip. The next active fault to the east is nearly purely extensional at another 3 mm yr −1 of extension. In southern Nevada a roughly east-west structure accommodates no more than 2 mm yr −1 of extension or strike slip (Fig. 13b). These low deformation rates are consistent with the Thatcher et al. (1999) hypothesis that little active deformation takes place within the Basin and Range.

Cascadia subduction and Pacific Northwest
For the full block model inversion we estimate spatially variable slip deficit rates on triangular dislocation elements representing the Cascadia subduction zone (McCrory et al. 2009) in addition to block rotation rates and fault slip rates (Fig. 15). In all models considered, we use a single smoothing parameter β = 1000. All three of the models considered fit the regional GPS observations with MRV <1.7 mm yr −1 (MRV = 1.6 mm yr −1 in the reference model). We estimate slip deficit rates of up to 22.7 mm yr −1 on the subduction zone (Fig. 15b). This increases to 28.2 mm yr −1 in the λ = 1150 model (Fig. 15c). This range is consistent with the coupling rates estimated by McCaffrey et al. (2007). Without modifying the smoothing parameter, the strain accumulation distribution on the subduction zone changes with the block geometry alone, and coupling rates on the subduction zone are particularly sensitive to the block geometry in the northwest Basin and Range. For example, a 21 mm yr −1 normal-sense patch on the southern tip of the subduction zone mesh in the λ = 1150 model (Fig. 15c) is likely related to fewer active structures in the Basin and Range in this model. All slip distributions estimate the highest slip deficit rates near the trench of the subduction zone, where slip is localized in two patches centred offshore northern California at 47 • N and offshore of the Olympic Peninsula at 41 • N, similar to previous estimates of strain accumulation on the subduction zone (e.g. McCaffrey et al. 2007McCaffrey 2009). Strain accumulation rates decay to half of the plate rate at ∼30 km depth. The slip distribution in our reference model is equivalent to a moment accumulation rate of a M W = 9.2 earthquake every 600 yr, consistent with palaeoseismic and recurrence interval estimates (Goldfinger et al. 2003;Satake et al. 2003) and previous geodetic imaging within block models (McCaffrey et al. 2007;McCaffrey et al. 2013).
In the reference model (λ = 850) and the λ = 650 comparison model, we identify five actively rotating blocks in the Pacific Northwest region: the Juan de Fuca plate, a coast range block, a large oblong block spanning Vancouver Island, eastern Washington, and northern Idaho, an eastern Oregon block and a southern Oregon/northeast California/northwest Nevada block (Figs 6 and  14). At the southern edge of the subduction zone and east of the Mendocino triple-junction, four roughly sub-parallel blocks accommodate northern California deformation. In Oregon and Washington, the coast range block is the only major block between the trench and northern Basin and Range. This block is bounded on the east by the Saint Helens Seismic Zone, the West Rainier Seismic Zone and the Straits of Juan de Fuca fault (McCaffrey et al. 2007). The northern Basin and Range is bounded on the east by the Bitterroot and Mission faults in western Montana, and to the south by several faults along the eastern Oregon-Washington border. The eastern Oregon block is separated from the interior Basin and Range block by a string of faults in western Idaho and the Brothers fault zone in southern Oregon. The block geometry is simpler in the λ = 1150 model, with a single coast range-northern Basin and Range block, a single Basin and Range block, and three blocks in northern California (Fig. 15c). As mentioned above, this simpler block geometry is associated with higher slip rates on the Cascadia subduction zone.
In all of the models considered in the Pacific Northwest, slip rates on block-bounding faults are low (0-5 mm yr −1 ; McCaffrey et al. 2007). Although we do not include faults in western Oregon that may accommodate additional right lateral strike-slip (e.g. McCaffrey et al. 2007), the lack of high fault slip rates especially in southern Oregon northern California provide negligible indications of additional deformation in this region. However, we do estimate right lateral strike-slip strain accumulation rates of 20 mm yr −1 on the subduction zone just north of the Mendocino triple-junction offshore northern California and southern Oregon (Fig. S2). These high right-lateral strain accumulation rates on the subduction zone may be consistent with deformation accommodated by instead block rotations in southern Oregon (McCaffrey et al. 2007), but which we have not included in these models.
Misfit in the Pacific Northwest remains consistently low (MRV < 2.5 mm yr −1 ) up to λ = 3150 (Fig. 5), suggesting that spatially variable slip on the subduction zone accommodates slip that is not otherwise accommodated on faults, and illustrates a trade-off between surface fault geometry and subduction zone slip. Therefore, the estimate of an initial slip distribution (Fig. 4) may contribute to final block geometries identified with the TVR algorithm. The influence of the initial slip estimate on block geometry is limited by the presence of a single coastal block over most of the subduction zone. However, active faulting in coastal Oregon and Washington not included in the initial dense block model may still map into spatially variable slip estimates on the Cascadia subduction zone.

DISCUSSION
The application of a total variation regularization (TVR) algorithm to an initially dense block model geometry allows us to identify the simplest block models that fit geodetic observations at a given resolution level. In this way, we have generated a suite of full plate boundary scale block models of the western United States by algorithmically identifying the tectonically most important block boundaries. At a residual level of 1.7 mm yr −1 , plate boundary deformation can be described with 30 blocks (λ = 850) in which deformation is concentrated along the San Andreas fault, San Jacinto fault, and the eastern California shear zone in southern California; the San Andreas and Calaveras systems in the San Francisco Bay Area; the Juan de Fuca plate boundary and Cascadia subduction zone in the Pacific Northwest and on discrete but low slip rate structures in the Basin and Range, as far east as the Wasatch fault.
The λ = 1150 model is the most quantized model that fits the horizontal GPS velocities with a MRV less than 2 mm yr −1 (Fig. 5), providing a bound on plate boundary complexity at this level of resolution. To fit observations better than 2 mm yr −1 , at > 4 mm yr −1 of slip is still required west of the San Andreas fault in southern California. Furthermore, western United States plate boundary deformation cannot be described at this resolution without an active Wasatch fault. The Wasatch fault remains included in the block model selection up to an MRV = 2.4 mm yr −1 (λ = 2100).
The reference 30-block model, can be decomposed into velocity components due to block rotations (Fig. 16a), and velocities due to elastic strain accumulation on locked faults and the Cascadia subduction zone (Fig. 16b), and illustrates the relationship between block rotations and the final interseismic velocity field. While the final velocity field (Fig. 16b) retains large-scale block features, especially along the central San Andreas fault and in the Basin and Range, sharp block boundaries in the denser fault systems of southern California and the San Francisco Bay Area are obscured by deformation due to elastic strain accumulation. These earthquake cycle elastic effects explain how the apparently smooth deformation field across much of the western United States is well described by slip on only a subset active of possible structures.
Also apparent in this decomposition are block rotations in the northeastern Basin and Range block. In general, we do not observe a systematic discrepancy between the geological slip rate constraints and our jointly estimated slip rates (e.g. Meade et al. 2013). Of the 32 geological estimates included in the λ = 850 reference model (on faults that remain active after TVR selection), 21 of our estimated slip rates agree within the geological uncertainty (Fig. 17). We estimate slip rates lower than the geological rate on 4 faults, and greater than the geological rate on 7 faults.
The eastern California shear zone is often cited as a region in which geological slip rates are systematically less than geodetically estimated rates (e.g. Dolan et al. 2007). At the geological study locations (Oskin et al. 2007;Oskin et al. 2008), we estimate slip rates of 0.6-2.0 mm yr −1 , comparable to the reported 724 E.L. Evans, J.P. Loveless and B.J. Meade Table S1; for details on how we have applied uncertainties, see Section 3). Selected fault labels shown in red. geological rates of 1-1.8 mm yr −1 . However, faults within 40 km of the Pisgah-Bullion fault (geological rate of 1 mm yr −1 , estimated rate of 0.6 mm yr −1 ), reach right-lateral rates of up to 11 mm yr −1 and left lateral rates of 15 mm yr −1 . We interpret these short-wavelength features as due to overfitting the combined GPS and geological observations, and not representative of present-day tectonics in the eastern California shear zone.
We estimate similar short-wavelength slip rate variation in the San Francisco Bay Area. This is perhaps apparent on the San Francisco peninsula segment of the San Andreas fault. The geological slip rate on the peninsula segment is 17 ± 4 mm yr −1 , and we estimate a slip rate of 23 mm yr −1 . This slip rate is partially accommodated by unrealistic left-lateral slip on the Zayante and San Gregorio faults to the west. This suggests that the inland Bay Area faults may not be consistent with slip rates required to allow high slip on the peninsula segment and maintain kinematic consistency across the boundary. Similarly, we estimate left lateral slip rates on the Valley Margin fault in the East Bay adjacent to the 23 mm yr −1 slip rate on the creeping segment of the San Andreas fault. The San Francisco Bay Area and the eastern California shear zone are both notable for hosting several closely (<20 km) spaced strike slip faults. In these regions, the overlapping interseismic strain accumulation signals due to locking on each of these faults inevitably leads to ambiguity and trade-offs in identifying fault slip rates (e.g. . At the GPS resolution and fault system density levels considered in this study, incorporating geological slip rate constraints does not necessarily resolve these ambiguities. Furthermore, observation density appears to correspond with fault system density in the reference model. In particular, both southern California and the San Francisco Bay Area feature the densest fault systems within the reference and comparison models, which may be due to large number of GPS stations and geological slip rate constraints in both southern California (415 stations and 28 geological rates) and the Bay Area (184 stations and 7 geological rates). To mitigate this effect, we may wish to additionally reweight the observations by the inverse of their spatial density, and/or by interpreting GPS-only block models in addition to joint solutions like those performed here. Additional improvements in the TVR block modelling methodology could be made by incorporating more complex tectonic and earthquake cycle processes by allowing spatially variable creep on select faults (e.g. ) and by accounting for time-dependent viscoelastic deformation across the boundary (e.g. Nur & Mavko 1974;Savage & Prescott 1978;Hammond et al. 2011;Hearn et al. 2013;Johnson 2013).

CONCLUSION
The kinematic block models of the western United States developed here are based on the idea that fault slip results from the differential motion of tectonic plates. We have treated the identification of the most active of these tectonic microplates as a simultaneous part of a total variation regularization optimization strategy that maximizes the goodness-of-fit to both interseismic GPS velocities and geological slip rates, while minimizing the number of unique microplate rotation vectors. This methodology allows construction of initially spatially dense block models and then effectively answers the question: 'How many kinematically distinct microplates (i.e. having a unique Euler pole) are required to explain the kinematic observables at a particular MRV?' For example, a three-block model (Pacific, North America and Juan de Fuca plates) with boundaries derived from estimation with a large differential rotation vector penalty, λ, fits the observed GPS velocities poorly with an MRV of 4.1 mm yr −1 , and localizes slip along the San Andreas fault. A reference model with 30 rotating blocks fits the GPS velocities with a MRV of 1.7 mm yr −1 . This approach to block model selection may be improved with the inclusion of additional deformation processes (e.g. spatially variable fault creep, time-dependent viscoelastic deformation), and serves as a step towards more algorithmic analyses of both plate boundary zone kinematics and earthquake cycle deformation.

A C K N O W L E D G E M E N T S
We thank editor Saskia Goes and two anonymous reviewers for their thoughtful reviews that greatly improved the content and clarity of the manuscript. Computational resources were provided and maintained by Harvard FAS Research Computing. This research was supported by funding from Harvard University and the Southern California Earthquake Center (SCEC), proposal number 12131. SCEC is funded by National Science Foundation Cooperative Agreement EAR-0106924 and U.S. Geological Survey Cooperative Agreement 02HQAG0008.

S U P P O RT I N G I N F O R M AT I O N
Additional Supporting Information may be found in the online version of this paper:   Table S1. Geological rates used in the TVR estimation and final block model solution.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.