Journal of Geophysical Research : Solid Earth A Network Inversion Filter combining GNSS and InSAR for tectonic slip modeling

Studies of the earthquake cycle benefit from long-term time-dependent slip modeling, as it can be a powerful means to improve our understanding on the interaction of earthquake cycle processes such as interseismic, coseismic, post seismic, and aseismic slip. Observations from Interferometric Synthetic Aperture Radar (InSAR) allow us to model slip at depth with a higher spatial resolution than when using Global Navigation Satellite Systems (GNSS) alone. While the temporal resolution of InSAR has typically been limited, the recent fleet of SAR satellites including Sentinel-1, COSMO-SkyMED, and RADARSAT-2 permits the use of InSAR for time-dependent slip modeling at intervals of a few days when combined. With the vast amount of SAR data available, simultaneous data inversion of all epochs becomes challenging. Here we expanded the original network inversion filter to include InSAR observations of surface displacements in addition to GNSS. In the Network Inversion Filter (NIF) framework, geodetic observations are limited to those of a given epoch, with a stochastic model describing slip evolution over time. The combination of the Kalman forward filtering and backward smoothing allows all geodetic observations to constrain the complete observation period. Combining GNSS and InSAR allows modeling of time-dependent slip at unprecedented spatial resolution. We validate the approach with a simulation of the 2006 Guerrero slow slip event. We highlight the importance of including InSAR covariance information and demonstrate that InSAR provides an additional constraint on the spatial extent of the slow slip.


Introduction
For a better understanding of what causes and triggers earthquakes, it is important to study all processes of the earthquake cycle. One aspect of this is the study of the spatial and temporal interrelation between the interseismic period, coseismic events and related post seismic signals, as well as aseismic slip processes such as slow slip events, which also change the surrounding stress field.
In the last few decades geodetic observations have proliferated with the development of dense permanent Global Navigation Satellite Systems (GNSS) networks such as the GPS Earth Observation NETwork in Japan, the Southern California Integrated GPS Network, the Pacific Northwest Geodetic Array, and the Sumatran GPS Array, and the acquisition of interferometric synthetic aperture radar (InSAR) data from a large variety of satellites. Multiple studies have used the high temporal resolution of continuous GNSS stations to model time-dependent processes including those of slow slip events [e.g., Cervelli et al., 2002;Segall et al., 2006;Schmidt and Gao, 2010;Radiguet et al., 2011;Bartlow et al., 2014;Ozawa et al., 2007;McGuire and Segall, 2003], post seismic slip [e.g., Hsu et al., 2006;Kositsky and Avouac, 2010;Miyazaki et al., 2004;Bedford et al., 2013], and transient deformation [e.g., Mavrommatis et al., 2014]. However, the spatial resolution is dependent on the local GNSS network and thus GNSS station distribution. In contrast, Interferometric Synthetic Aperture Radar (InSAR) has a much finer spatial resolution, on the order of meters, but is limited to longer time scales, with acquisitions every few days at best, and is only sensitive to deformation in the direction of the radar line of sight. Because of complementary advantages, GNSS and InSAR are often used in a joint framework [e.g., Pritchard et al., 2002;Simons et al., 2002;Wright et al., 2004].
One way to retrieve the time-dependent history of fault slip is to invert all observations at all epochs simultaneously for slip at depth. This can become data and memory intensive, especially when considering vast amounts of continuous GNSS and InSAR data, the latter of which can be a few millions of observations for a single track and epoch. Strategies exist to decrease the amount of InSAR data but might not be sufficient, In our study we focus on the extension of the NIF (version by Bartlow et al. [2014], which is an expansion of Segall and Matthews [1997]) to combine GNSS and InSAR observations. The NIF implementation uses Kalman forward modeling and a backward smoothing operation. After completing both operations, all geodetic observations will provide a constraint on the slip estimates at all epochs. This is different to the GA-LKF approach [e.g., Shirzaei and Walter, 2010] that runs forward in time and is iterated until convergence is reached. The strength of the NIF thus lies in the complementary InSAR and GNSS data sets, where GNSS provides a high temporal resolution and InSAR gives high spatial resolution. We describe and implement the methodology required to include InSAR in the NIF, combined with GNSS. We then demonstrate the procedure on a synthetic simulation of the 2006 Guerrero slow slip event.

Network Inversion Filter
To model time-dependent fault slip, the relationship between fault slip, s, and geodetic observations of surface displacements, d, is combined with a stochastic description of how slip evolves over time. The relationship between slip and the displacement at the surface follows from elastostatic Green's functions [e.g., Okada, 1985;Thomas, 1993]. However, surface displacements as observed by geodetic techniques will be contaminated by processes including nontectonic deformation, such as motions introduced by soil compaction. In addition, GNSS reference frame corrections, InSAR orbit errors, and atmospheric delays introduce additional apparent deformation that must be accounted for.

10.1002/2015JB012638
The definition of the transition matrix depends on the stochastic model. The network inversion filter (NIF) as proposed by Segall and Matthews [1997] is designed to detect the departure of slip from steady state, which we define to be the interseismic rate. The prior assumption is that the slip acceleration is close to zero, and can be modeled as a white noise process, with scale parameter as ∼ ( 0, 2 ) . This allows cumulative slip s since time t 0 to be written as where v is the interseismic slip rate and W t the accumulated slip deviated from the interseismic rate. The accumulated transient slip is the integral of a random walk processẆ or twice the integral of a white noise process with variance 2 . The scale parameter constrains the temporal smoothing of the slip.Ẇ is the transient slip rate. For a more complete description of the theoretical basis for the NIF see Segall and Matthews [1997].
Below, we elaborate on the observation equations for GNSS (section 2.1) and InSAR (section 2.2) and describe how the state variables are assumed to change in time. Finally (section 2.3), we combine both GNSS and InSAR and derive the full observation matrix H k , the observation variance-covariance matrix R k , the state transition matrix T k+1 , and the process noise variance-covariance matrix Q k+1 . These four matrices completely define the linear system.

GNSS Observation Equation
The GNSS surface displacements, d GNSS , can be written as [e.g., Segall and Matthews, 1997] where x describes the station location, t GNSS 0 is the start of the GNSS time series, G are Green's coefficients relating slip s to surface displacements, s t GNSS 0 the slip since t GNSS 0 , and  are the local GNSS benchmark motions for each component (East, North, and Up) and every station, modeled as a Brownian random walk with scale as  = ∫ t 0 dw [Langbein and Johnson, 1997]. Note that the benchmark motions represent spatially incoherent GNSS network displacements and therefore should not include displacements related to Gs. F is the GNSS reference frame error, where F is a linearized Helmert transformation [e.g., Miyazaki et al., 2003;Mavrommatis et al., 2014].  is a vector containing the coefficients of the Helmert transformation (translation, rotation, and scale factors for each component). We assume F to be an identity matrix, and let the 2 control the variance of the Helmert transformation coefficients. GNSS are the GNSS observation errors ∼ (0, GNSS (t)), where GNSS (t) is the GNSS observation variance-covariance matrix at time t.

InSAR Observation Equation
The InSAR surface displacements, Δd InSAR , are the difference in the radar line-of-sight displacements between two acquisition times, t InSAR 0 and t, and for which the observations are with respect to an arbitrary reference area or pixel [e.g., Hooper et al., 2012;Bekaert et al., 2015a] where t InSAR 0 refers to the acquisition time of the first SAR image and P = [ ] ⊤ is a planar correction to account for the long wavelength orbit errors in the interferogram between t InSAR 0 and t; x 1 and x 2 are the InSAR observation geocoordinates converted in a local reference frame (x 1 , x 2 in unit of meters) where the local origin is defined with respect to the center of the InSAR study area. The variation of the planar coefficients in time is assumed to be uncorrelated between interferograms. InSAR are the InSAR observation , where InSAR (t) is the interferogram observation variance-covariance matrix at time t of the noise. Note that this is the variance-covariance matrix in space and not in time.

Joint GNSS and InSAR Observation and State Transition Equation
For the derivation of the full observation matrix H k , assume that the number of GNSS stations, N s , does not vary in time. Likewise for the number of InSAR pixels, N p , which is consistent over N I interferograms. Note that consistent observations over time are not required in order to execute the NIF. After the general derivation of the observation matrix we elaborate how to modify H k when certain GNSS observations and InSAR pixels are not available for a given epoch. In our derivation, we also limit the description below to one GNSS network and InSAR track but expand H k to include multiple InSAR tracks in Appendix A. Combining the GNSS observations, equation (4), with those of InSAR, equation (5), in the observation equation (1), including Laplacian smoothing, at time t k results in where Δd InSAR is a ( N p × 1 ) vector of InSAR line-of-sight displacements and d GNSS a ( 3N s × 1 ) vector containing the three-component GNSS displacements since t 0 . Here we choose to minimize the Laplacian of the interseismic slip rate (v), transient slip (W), and transient slip rate (Ẇ) separately. Alternative options can be included, such the slip (s = vt + W) and/or slip rate (̇s = v +Ẇ). Assuming that the fault is modeled using N d dislocation patches, the state vector comprises an interseismic rate vector v, an integrated random walk vector W (cumulative slip deviation from the interseismic), and a random walk vectorẆ (transient slip rate)-all of length N d . Additional terms are defined in equations (4) and (5). The length of the state vector does not change and is identical in the state observation and state transition equations. As the transient slip rate does not have an influence on the observations, the column in the observation matrix consists of zeros, except for the Laplacian smoothing included through the pseudo observations. We assume that the InSAR network is defined or inverted with respect to the first acquisition at time t InSAR 0 . By doing so, we are able to estimate the reference slip with respect to t InSAR 0 . In cases of isolated networks in time on a single track, the subnetworks can be regarded as "new" tracks (see Appendix A on how to include multiple InSAR data sets). Initially, the reference slip is assumed to be zero or a prior slip is assumed, which is updated to s InSAR We include the reference slip at the time of the master acquisition by including a set of pseudo observations in the observation equation as 0 = vt k + W k − s InSAR 0 . This is not shown in equation (6) and only applies at t k = t InSAR 0 . Note that the initialization of the state vector constrains the reference slip in case no GNSS observations are present prior to the master acquisition. At epochs where the GNSS and InSAR data are not available, the observation equation (2.3) is modified by deleting those rows of the data vector and observation matrix corresponding to the missing GNSS stations or InSAR locations.
The observation variance-covariance matrix follows from the combination of GNSS, InSAR, and the weight of the Laplacian smoothing, 2 i , assuming all to be uncorrelated with each other as where 1 and 2 are scale parameters to account for relative scaling between GNSS and InSAR and to account for errors in the model. We assume that the observation variance-covariance matrices are well described and do not require a relative weighting, which allows us to simplify equation (7) with = 1 = 2 . Note that the hyperparameter still allows us to account for model errors. For simplicity, we assume the weight of the Laplacian smoothing to be the same for the interseismic rate, slip, and slip rate; thus , = 1 = 2 = 3 . When GNSS or InSAR data are not available for a given epoch, the data covariance matrix, equation (7), is modified by deleting both columns and rows corresponding to those missing GNSS stations or InSAR locations.
In the state transition equation (2), the interseismic slip rate is by definition constant in time and does not, therefore, change from epoch t k to t k+1 . The transient slip (integrated random walk) at the new epoch, W k+1 , follows from that of the previous epoch combined with the integration of the random walk between both epochs as W k+1 = W k + ( t k+1 − t k )Ẇ k . As indicated before, the GNSS benchmark motion follows a random walk. The GNSS reference frame correction is assumed to be independent from epoch to epoch [e.g., Miyazaki et al., 2003] , where is a scale parameter of the reference frame 10.1002/2015JB012638 variance-covariance matrix I  . Similarly, the InSAR orbit (plane) is assumed to be white noise ∼ ( 0, 2 I  ) . In our approach we invert all interferograms to a common master. While interferograms are correlated in time, we do not account for this in our model. The full state transition equation can be written as with 0 a zero matrix, I an identity matrix, both of size For the process noise variance-covariance matrix, Q t k , we follow, e.g., Segall and Matthews [1997]; Segall et al. [2000] for the interseismic slip rate, integrated random walk, random walk, and GNSS benchmark motion; and Miyazaki et al. [2003] for the GNSS reference frame errors. We apply the same methodology for the InSAR orbit (plane) with error distribution ∼ ( 0, 2 I  ) . Like the interseismic slip rate, the reference slip does not change in time. The full process noise variance-covariance matrix Q t k is therefore . Like before, 0 and I are the zero and identity matrices. The size of I  is (3 × 3), while for I  the size depends on the parameters included; for example, in case of translation, rotation, and scaling, it has a size of (7 × 7).

Kalman Filter and Backward Smoothing Procedure
The full description of the Kalman forward filtering and backward smoothing approach is contained in Segall and Matthews [1997]. During the forward filtering step, the stochastic model is used to predict the state at the next epoch,X k+1|k , based on the state of the current epoch,X k|k . We use the notationX k+1|k , which readsX at k + 1 given k and is the state at epoch k + 1 given the data up to epoch k. The prediction is described by the state transition equation (2) This is followed by an update to the estimated state, by including the observations of epoch k + 1 aŝ The notationX k+1|k+1 now readsX at epoch k + 1 given the data up to epoch k + 1. K k+1 is the Kalman gain at epoch t k+1 defined as The prediction and update procedure is applied in sequence to the whole time series until the observations of the last epoch, N, are updated, yieldingX N|N . At this stage, the estimated state is conditional on the data up to that epoch. After performing the backward smoothing operation, where the same recursive prediction and update structure is used as during the forward Kalman filtering but with time reversed (i.e., starting from the last epoch), all geodetic observations are used in constraining the state vector at all epochs, yieldingX k|N for all k [Rauch et al., 1965;Segall and Matthews, 1997].
To initiate the Kalman filter, an a priori estimate of state vector given no data,X 1|0 , is assumed, which describes the state at the initial epoch without any constraint from the data. The corresponding uncertainties are specified in 1|0 . Larger uncertainty can be attributed when good a priori knowledge of the state vector is lacking.
The NIF requires the specification of the hyperparameters for spatial and temporal smoothing. These can be selected through a maximum likelihood grid search, as proposed in the original NIF [Segall and Matthews, 1997], and applied in later studies [e.g., Segall et al., 2000;Bartlow et al., 2011Bartlow et al., , 2014. Unlike previous studies, we include both observations and pseudo observations in the computation of the maximum likelihood. Including the pseudo observations introduces a trade-off between fitting the observations and smoothing, which otherwise would not be included, leading to a higher likelihood for a rougher solution. Excluding pseudo observations would therefore lead to a rougher maximum likelihood solution. Once the hyperparameters have been selected, these can be ingested in a new NIF run.

Comparing the NIF With the Earlier Implementation
We modified the version of the NIF by Bartlow et al. [2014], which works with GNSS data only, to include also InSAR data. While the main changes are related to the InSAR data component, we also include the interseismic rate directly in the inverse problem, whereas Bartlow et al. [2014] correct the GNSS observations for the interseismic slip rate prior to ingestion in the NIF. To assess the impact of our modifications, we compared the results of our NIF version and that of Bartlow et al. [2014] when inverting GNSS observations of the 2011 Cascadia Slow Slip Event (SSE). This data set is distributed as test data for the earlier implementation of the NIF. Overall, we find that our implementation does not significantly alter the results obtained from the previous implementation; average transient slip differences are <1 cm, and slip rate differences are negligible (∼0 cm/d).

Synthetic Simulation of the Guerrero SSE
Our test data set reflects the tectonic setting in southern Mexico ( (1 year duration). We simulate the SSE, guided by the time-dependent GNSS modeling results of Radiguet et al. [2011], who solved for the slow slip source time function. We adopt a similar source time parameterization scheme [Liu et al., 2006], where the location of slow slip initialization and the 0.8 km/d slip propagation velocity is from the best fit GNSS model by Radiguet et al. [2011]. For simplicity we fix the slow slip rise time, which is estimated to vary between 160 and 200 days [Radiguet et al., 2011], to be 183 days. For the maximum accumulated slip over the duration of the 2006 SSE, we use the results from Bekaert et al. [2015a], which are similar to other studies of the same event [e.g., Vergnolle et al., 2010;Radiguet et al., 2011;Cavalié et al., 2013]. A full mathematical description of our time-dependent slow slip model as a function of the local slow slip start time, local rise time, and local cumulative slip magnitude is contained in Appendix B. Spatial maps of the local start time, rise time, and cumulative slip are shown in Figure S5 in the Supporting Information. We define the interseismic loading using a back slip formulation [Savage, 1983], where the slip deficit follows from the multiplication of the time with the MORVEL (Mid-Ocean Ridge VELocity) plate convergence rate of 6.1 cm/yr [DeMets et al., 2010] and a simulated inter-SSE coupling. The latter only varies with depth; we assume that the fault is freely slipping at shallow depths (coupling = 0), fully locked in the seismogenic zone between ∼15 and 25 km (coupling = 1), and that there is a smooth transition region from 25 to 50 km depth, below which the fault is freely slipping.
We use the time-dependent cumulative simulated slip (interseismic and slow slip) on the subduction interface in combination with Green's functions from triangular dislocations [Thomas, 1993] to infer the displacements over time at the surface. We use the real GNSS station distribution that was installed in 2006 but limit ourselves to seven stations, with six in the far field (magenta markers in Figure 1). We do this to show the impact of including InSAR more clearly in a typical region with sparse GNSS coverage. GNSS observations are simulated every 3 days, with random white noise of 1 mm uncertainty for the horizontal displacements and 2 mm for the vertical displacements. We assume no additional local GNSS station motion or reference frame motion but will  [Pardo and Suárez, 1995;Melgar and Pérez-Campos, 2011;Pérez-Campos and Clayton, 2014].
still solve for this in our NIF. A time series with the north, east, and up components is shown for selected GNSS stations in Figure 2a. The location of our simulated InSAR track corresponds to the location of the descending Envisat track 255 (red polygon in Figure 1) used in earlier slow slip studies over the region [e.g., Bekaert et al., 2015a;Cavalié et al., 2013]. We fix the master SAR acquisition to be on 1 June 2005, with a simulated SAR acquisition every 2 months. We generate 10 interferograms, and include orbit errors by simulating the addition of a bilinear plane, for which the coefficients are summarized in Table 1. We incorporate the effects of   Lohman and Simons [2005], with L c the exponential range and rn the standard deviation of the uncorrelated noise. residual atmosphere delays by including spatially correlated noise according to Lohman and Simons [2005], which varies in magnitude and correlation length for each interferogram as summarized in Table 1. Figure 2b shows the individually simulated components of the interferogram and the interferograms as input to the NIF. Our simulated interferograms have similar signal magnitude as processed Envisat data over the region [Bekaert et al., 2015a].

Results
We compare the NIF results when inverting for the GNSS observations only and when jointly inverting GNSS and InSAR observations. InSAR covariance is often neglected in slip inversion studies. This should be avoided as spatially correlated atmospheric noise will be treated as if it were a signal. The importance of covariance information has also been highlighted in other studies [e.g., Lohman and Simons, 2005;Hetland et al., 2012;Bekaert et al., 2015a;Agram and Simons, 2015;Fattahi and Amelung, 2015]. To demonstrate the impact of InSAR covariance on the estimated slip, we also include a comparison between results when the InSAR covariance is included in the inversion and when it is neglected. Our variance-covariance matrix is based on an exponential covariance function [Wackernagel, 2003] with the range and variance the same as those that were used when simulating the spatially correlated atmospheric noise (Table 1). We do not change the initialization of the state vector parameters between the different cases. At t = 0, we assume the interseismic slip rate to be 0 cm/yr  Figure 4b, the cumulative slip in Figure 5b, and the slip rate history in Figure 6. The red circle marker in Figures 3b and 3c corresponds to the maximum likelihood solution when using only GNSS (Figure 3a). with a 6 cm/yr uncertainty, the accumulative slip to be 0 cm with a 1 cm uncertainty, and the transient slip rate to be 0 cm/d with an uncertainty of 10 −8 cm/d. These parameters are chosen assuming that the observations start in an inter-SSE period. In our modeling we account only for reference frame translation (east, north, and up components). A priori knowledge of GNSS reference frame errors [e.g., Wdowinski et al., 1997] and local station motion errors [e.g., Dmitrieva et al., 2015] can be used to set their uncertainty in the NIF. However, we choose the reference frame uncertainty (1 mm) and local station motion uncertainty (1 mm yr −0.5 ) by trial and error such that local variations are allowed over time, but no tectonic leakage can be observed in time. This approach is needed for two main reasons: (1) the SSE is observed at only one GNSS station and therefore looks like an incoherent motion, while in reality it is a tectonic signal and (2) all GNSS stations are located on one side of the fault (hanging wall), and therefore, the interseismic slip rate is quasi-uniform and will be partly mapped into the GNSS frame motion. We set the uncertainty of the parameters defining the orbital ramp to be high (e.g., 10 m and 10 m/km), which allows us to accurately fit the best bilinear plane between GNSS and InSAR. The spatial and temporal smoothing hyperparameters vary between the different cases and are selected though a maximum likelihood grid search [e.g., Segall and Matthews, 1997;Segall et al., 2000;Bartlow et al., 2011Bartlow et al., , 2014. As all GNSS observations are in the reference frame of the hanging wall, we solve for the interseismic loading in a back slip framework [Savage, 1983]. The transient slip has the opposite sense to the interseismic slip deficit. To avoid leakage of the interseismic slip rate into the transient slip rate and vice versa, we enforce a negativity constraint on the interseismic slip rate consistent with the back slip formulation while enforcing a positivity constraint on the transient slip rate. Identical to the approach by Bartlow et al. [2011], we use the PDCO (primal-dual interior method for convex objectives) package [Saunders, 2015] to optimize the state vector after imposing the sign constraint. This approach works in two steps. First, the interseismic slip not well distributed, with only the ACAP station clearly capturing the on their sign constraint [Simon and Simon, 2010]. Second, we use the PDCO package [Saunders, 2015] to optimize the state vector based on a covariance-weighted L2-norm minimization of the difference between the constrained and unconstrained rates. Due to computation intensity, we do not run the PDCO optimization for the state vector as part of our hyperparameter grid search. However, after the maximum likelihood estimation (MLE) of the hyperparameters, we compute the MLE slip solution including the PDCO optimization. Below, we report on the results from the maximum likelihood solution for the three defined cases. Figure 3a shows the results of the hyperparameter grid search, with a maximum likelihood value for the spatial smoothing parameter of 2 ∕ 2 = 24 ⋅ 10 −8 , temporal smoothing parameter of 2 ∕ 2 = 0.24, and estimated data variance-covariance scaling of̂= 0.97. As expected, the latter is close to 1, indicating that the inversion is estimating the model to within expected errors. The results of the maximum likelihood solution after forward filtering and backward smoothing are shown in Figure 4 for the inter-SSE locking rate, Figure 5 for the cumulative estimated slow slip, and Figures 6-8 for the transient slip rate history.

GNSS Only
The coastal GNSS stations are located approximately above the locked part of the subduction interface, Figure 4a, which is loaded in our simulation at a rate of ∼6 cm/yr. We are capable of retrieving this same peak magnitude of 6 cm/yr. We find that the locking rate on the fault is well estimated for patches close to (∼20 km) the coastal GNSS sites, with a residual that falls within the ∼1.5 cm/yr estimated uncertainty (Figure 4c). While a similar strike-parallel pattern can be observed, we find the peak distribution to be located ∼10 km shallower than that in the input simulation. We observe some smearing farther downdip, with an average residual of ∼2.5 cm/yr, which is likely due to the imposed smoothing, to compensate for the overestimation near the trench. The poor resolution is not surprising given the very sparse GNSS coverage.
Our estimated cumulative slip, Figure 5b, is a smeared version of the simulation, Figure 5a. This is to be expected, as the GNSS network is not well distributed, with only the ACAP station clearly capturing the SSE surface displacements. We find the peak slip to be underestimated by ∼7 cm, with an average estimated uncertainty of ∼2 cm (Figure 5c). We also find some updip slip residuals introduced by the misestimation of the inter-SSE loading. The latter also causes the slip uncertainty to grow over time.
We find the estimated transient slip rate (Figure 6, second row) to correlate in space and time with the simulation ( Figure 6, first row). Supporting Information Figure S6 shows the same but for a 9 day interval without averaging. As for the cumulative slip, we find the estimated transient slip rate to be smoother than that of the  Figure 9 shows the variance and covariance history between v, W, andẆ (red line) for the triangular dislocation with peak slow slip and for another away from the slow slip region. As expected, we find that the interseismic slip rate is negatively correlated with the cumulative transient slip and transient slip rate. Halfway through the observation period, the covariance between v andẆ flattens. This could mean that once v is well established by the pretransient data it does not trade off as much with transient slip. The correlation between v and W keeps increasing (in absolute value) because W integrates any uncertainty inẆ. The interseismic slip rate and its uncertainty are fixed based on the estimate from the last epoch and do not vary when performing the backward smoothing operation.  Original observations and those estimated are shown in the Supporting Information ( Figure S1). The modeled surface displacements fit the observations well, with a mean residual of around 0 mm and an accuracy similar to that of the simulated observations: ∼1 mm for horizontal and ∼2 mm for the vertical components. Figure 3b shows the results of the grid search for the hyperparameters. The estimated maximum likelihood value for temporal smoothing hyperparameter ( 2 ∕ 2 = 0.8) is double that of the GNSS case. Compared to GNSS case, a much rougher spatial solution is preferred ( 2 ∕ 2 = 5⋅10 −6 ). We find a lower value for the estimated data variance-covariance scaling,̂= 0.89, which implies that the model overfits the data.

GNSS and InSAR (Neglecting Covariance)
As InSAR covariance information is not included, each of the ∼10 3 InSAR observations is treated as if it is an independent observation. A much rougher solution is therefore found for the maximum likelihood solution. In reality much of this roughness is fitting the spatially correlated InSAR noise from the residual atmosphere. This results in apparent noise signals that propagate into the slip model. The estimated inter-SSE loading rate ( Figure 4b) shows large residuals underneath the InSAR track, especially toward the downdip extent of the subduction interface (up to 8 cm/yr).
The cumulative transient peak slip is overestimated by ∼3 cm, Figure 5b. Also, large errors relative to the true input slip with an average of ∼12 cm can be observed at other locations, especially downdip of the slow slip region. Because the solution is rougher than the GNSS-only solution, we also find larger uncertainties for the estimated cumulative slip. The smallest uncertainty of ∼4.5 cm can be found underneath the InSAR track and GNSS station locations, with larger values with an average of ∼8 cm elsewhere.
The slip rate history (Figure 6, third row) shows much more temporal variation over the whole time period, even when there is no SSE taking place. These fluctuations reach an average magnitude of ∼0.07 cm/d. Also, the corresponding uncertainties are larger than when inverting for GNSS alone (Figure 7, second row).
We find larger magnitudes for the variance and covariance history (blue line in Figure 9) than when using GNSS alone. We also observe a correlation between changes in the covariance history and those times when InSAR data were ingested into the NIF.
We find that the GNSS residuals are within the uncertainty bounds of the original simulated observations (mean residual around 0 mm); see Supporting Information Figure S2a. For the InSAR data, we find an average root-mean-square error (RMSE) of ∼3 rad (or ∼1.3 cm). As we simulated the original tectonic signal and the InSAR orbit error (phase ramp), a direct comparison can be made with the NIF estimates. Supporting Information Figure S2b gives the simulated tectonic signal and the misfit of the estimation. We find a mean RMSE of ∼0.57 rad (or ∼0.25 cm). The last three interferograms (October 2006, November 2006, and February 2007 show the largest residuals, with an average of ∼1 rad (or ∼0.45 cm). Over time, the spatial pattern of the residuals is enhanced, illustrating that misestimation is integrated over time. We find that the

GNSS and InSAR (Including Covariance)
By including the full InSAR variance-covariance information, we find that the combined GNSS and InSAR solution is more similar to the GNSS-only case. This is expected as the InSAR observations add more information about the spatial extent of the slow slip surface observations. However, the InSAR contribution is downweighted and provides only a lower level constraint [e.g., Bekaert et al., 2015a], as the simulated slow slip signal has a correlation length similar to that of the simulated residual atmosphere (i.e., ∼2-50 km). By including the covariance, we find that the estimated data variance-covariance scaling is close to 1,̂= 0.99. We also find similar hyperparameters to the GNSS-only case.
Overall, an improvement can be observed when jointly inverting GNSS and InSAR, compared to the GNSS-only case, while accounting properly for the InSAR covariance information. Apparent noise signals, introduced in the slip rate history when neglecting the covariance information, are now suppressed (Figure 6). We find a 10.1002/2015JB012638 significant improvement in the location and estimation of the peak cumulative transient slip ( Figure 5), misfit of <1 cm compared to ∼7 cm for the GNSS-only case, and with a similar uncertainty of ∼2 cm. The inter-SSE loading is slightly different to the GNSS case ( Figure 4). We find a residual rate of up to 4.5 cm/yr underneath the InSAR track, downdip of the slow slip region. This is likely introduced to compensate for the smeared downdip slow slip signal. When neglecting InSAR covariance, this residual was significantly larger (up to ∼8 cm/yr).
We find that the variance and covariance history (green line in Figure 9) tends more toward the GNSS case. However, a correlation can still be observed between the time of InSAR data ingestion and changes in the covariance history.
As before, we find small residuals for the GNSS and InSAR (Supporting Information Figure S3). GNSS residuals are within the uncertainty bounds of the simulated data. For the InSAR, we find the estimated tectonic signal and orbital errors with similar magnitudes as when neglecting the InSAR covariance. The misfit between the tectonic simulation and the estimation indicates underestimation of the slow slip signal.

Discussion
Our NIF results are a reasonable approximation of the simulated tectonic slow slip signal. We find that the InSAR observations provide an additional constraint on the slow slip signal at the fault interface, which is underestimated when inverting sparsely distributed GNSS observations only. After performing the backward smoothing, the NIF is capable of estimating the spatial variation of the transient slip rates at intervals shorter than the InSAR observations (see Supporting Information Figure S6). When comparing with our simulation, we find a similar transient slip rate nucleation zone (location difference <25 km). From visual inspection we find that our NIF results are capable of capturing the propagation of the slow slip over time. In general, we find the transient slip rates to be slightly smoother and consequently with a smaller peak magnitude than simulated. A more quantitative approach, which could be investigated in future studies, would be to parameterize the slip history for each of the dislocation patches and then to invert directly for the slow slip nucleation and propagation speed [e.g., Radiguet et al., 2011].
Our NIF implementation is limited to a single spatial smoothing hyperparameter. We believe that it could be more appropriate to include a separate smoothing hyperparameter for the inter-SSE locking rate and also smooth differently in the dip and strike direction. For example, at subduction zones the interseismic slip rate is expected to vary with depth with a shorter wavelength than in the along-strike direction. Including an additional hyperparameter in our current NIF, implementation comes at the cost of an extra dimension to the grid search. In addition, but beyond the scope of our work, is the inclusion of model uncertainties due to, for example, fault geometry and elastic structure. For more information on how to handle these in a Bayesian framework see, for example, Duputel et al. [2014].
With a smaller grid spacing, and with more hyperparameters to solve for, the grid search becomes a time-consuming operation. An alternative approach to the grid search is included in the Extended NIF by McGuire and , where the hyperparameters are included directly in the state vector. As in our study, the hyperparameters remain constant over time. Assuming that constant hyperparameters can be a limitation in a complex time series with a mixture of interseismic, coseismic, post seismic, aseismic, and slow slip processes. In particular, the rapid displacement change during a short-term SSE might be suppressed due to temporal smoothing. Developments by Fukuda et al. [2004Fukuda et al. [ , 2008 allow for a time-varying temporal smoothing parameter. This is achieved by including the temporal smoothing hyperparameter as stochastic variables using the Monte Carlo mixture Kalman filter or the hierarchical Bayesian state space approach. None of the above NIF modifications currently include InSAR capability. We believe that our presented methodology can be adopted straightforwardly into the other methods. Shirzaei and Burgmann [2013] used the approach of Shirzaei and Walter [2010] to model time-dependent creep using surface creep data and an InSAR time series. The latter was preprocessed to integrate both Envisat and ERS data time series together and also filtered in space and time to reduce the contribution of atmospheric noise. In their approach a genetic algorithm is used in combination with a linear Kalman filter. Whereas we use a backward smoothing step, Shirzaei and Walter [2010] iterate their Kalman filter several times by using the estimates at the last epoch as initialization. While no further details are specified on the implementation of their method, further differences are also likely, such as how the relative nature and orbits of InSAR are handled and how the signal covariance of the observations, if included, is handled in time.
In our simulation, significant spatially correlated atmospheric noise (simulated variance between 1 and 3 cm 2 , and correlation lengths up to 50 km) is added. InSAR covariance information is often neglected in tectonic slip inversion studies. This should be avoided, as spatially correlated noise will be treated as signal. This becomes of special importance in studies where the tectonic signals have a similar correlation length as the noise, as is the case for our simulation here. Recently, different efforts and promising progress have been made on the quantification and construction of InSAR variance-covariance matrices [e.g., González and Fernández, 2011;Agram and Simons, 2015;Fattahi and Amelung, 2015;Bekaert et al., 2015a]. Future modeling studies should aim to include InSAR covariance information as part of regular inversion procedures.
The quality of the InSAR observations is expected to improve with more regular SAR acquisitions. Our simulation included a SAR acquisition every 60 days. With Sentinel 1 acquisitions are currently every 12 days, which will further decrease to every 6 days once Sentinel 1B is launched. With a swath width of ∼250 km and a spatial resolution of 5 × 20 m, large areas can be covered at high resolution and with short repeat times. Having a larger amount of independent InSAR observations will lead to a further reduction in the influence of atmospheric noise. In addition, atmospheric InSAR noise can be further reduced using tropospheric correction methods [e.g., Doin et al., 2009;Jolivet et al., 2011;Bekaert et al., 2015bBekaert et al., , 2015cFattahi and Amelung, 2015]. Being able to ingest InSAR data for time-dependent modeling in the NIF could provide invaluable information in a variety of applications, such as volcano inflation and deflation, and fault creep, as well as coseismic and post seismic events.

Conclusions
Vast amounts of geodetic data have been acquired over the last two decades. The acquisition rate will continue to increase exponentially with further expansion of GNSS networks and the acquisition of InSAR data from Sentinel 1, NISAR, and ALOS 2. Simultaneous inversion of all geodetic data to solve for the time-dependent history of fault slip is a computationally intensive process, for which an alternative is the Network Inversion Filters. Different versions exist of this method, but to date, none of these include Interferometric Synthetic Aperture (InSAR) observations. In this study, we have provided and applied the network inversion filter methodology to include InSAR. To validate the approach, we simulated the 2006 Guerrero SSE for a subset of the existing GNSS sites (mainly far field) and for the descending Envisat track 255. We find that GNSS can retrieve the cumulative SSE at approximately the same location as the input simulation, but with a smaller peak slip, due to spatial smoothing. Inclusion of high-resolution InSAR further improves the recovery of the transient slip. InSAR covariance is often neglected in slip inversion studies. This should be avoided, as spatially correlated atmospheric noise will be treated as if it was a signal, introducing apparent slip signals at depth. We compare the joint GNSS and InSAR inversion while neglecting and including InSAR covariance information. When including InSAR covariance, we find a solution that has similar smoothing hyperparameters to that of GNSS. We find that InSAR provides an improved constraint with its high-resolution observations above the slow slip region. Our study demonstrates the use of InSAR data to retrieve time-dependent slip, which can provide invaluable information for a wide variety of applications.  normal distribution; random walk-scale parameter controlling temporal smoothing; scaling parameter of the data observation variance-covariance matrix; scale parameter of the Brownian random walk of the local GNSS benchmark motion; scale parameter of the GNSS reference frame variance-covariance matrix; N s , N p , N I number of GNSS stations, InSAR pixels on a single track, and InSAR interferograms on a single track t GNSS 0 , t INSAR 0 GNSS and InSAR reference time; t time; v interseismic slip rate; W integrated random walk or cumulative slip deviated from the steady state interseismic slip rate;