Oblique slip on the Puysegur subduction interface in the 2009 July M W 7.8 Dusky Sound earthquake from GPS and InSAR observations: implications for the tectonics of southwestern New Zealand

The M W 7.8 Dusky Sound earthquake of 2009 July 15 was the largest earthquake in New Zealand in the past ∼ 80 yr and is the only major subduction interface earthquake in the New Zealand historical record. We have estimated the coseismic and early post-seismic slip distribution in the earthquake by inversion of GPS and differential interferometric synthetic aperture radar (DInSAR) observations. We show that slip during the earthquake was highly oblique and essentially in the Australia–Paciﬁc relative plate motion direction. Failure occurred on a ca . 80 by 50 km patch of the subduction interface with maximum slip 5–6 m. In this region, the offshore Alpine Fault carries up to ∼ 30 mm yr– 1 of slip before terminating at Resolution Ridge near the southern end of the 2009 slip patch. The southernmost part of the Alpine Fault trace lies near the up-dip end of the 2009 slip patch and the fault probably exists only within the crust of the over-riding plate in this region, terminating at shallow depth where it meets the subduction interface. There is no evidence for rupture of the Alpine Fault during the 2009 event. However, Coulomb stress changes from the subduction earthquake have brought this section of the Alpine Fault closer to failure. In contrast to many subduction zones undergoing oblique subduction, where the slip in major earthquakes is partitioned into more-or-less pure thrust earthquakes on the subduction interface and strike-slip earthquakes in the backarc, this is a case where the majority of the subduction interface slip is not partitioned, with only a shallow convergence zone west of the strike-slip Alpine Fault undergoing contraction approximately normal to subduction zone strike. North of the 2009 earthquake region the Alpine Fault lies further east relative to the subduction interface and partitioning between strike-slip on the Alpine Fault and trench-normal thrusting on the subduction interface takes place in a more typical fashion. Our results provide the ﬁrst clear demonstration of non-partitioned oblique slip between the Australian and Paciﬁc plates offshore of southern New Zealand.


I N T RO D U C T I O N
The 2009 July 15 M W 7.8 Dusky Sound earthquake ruptured a ∼80 km long section of the Puysegur subduction interface, where the Australian Plate is obliquely subducted beneath southwestern New Zealand ( Fig. 1; Fry et al. 2010). It was the largest earthquake in New Zealand since the 1931 February 3, M W 7.8 Napier earthquake and is the largest event ever recorded at the Puysegur subduction zone.
Our understanding of the Puysegur subduction zone has developed substantially since early workers identified its Neogene age  ; Barnes et al. 2002); the offshore Alpine Fault (Barnes et al. 2005) and the onshore Alpine Fault (Sutherland et al. 2006). The white arrow shows the present-day relative plate motion from Beavan et al. (2002); the MORVEL (DeMets et al. 2010) plate motion model has the same rate with a direction 4 • closer to north. Resolution Ridge is the southernmost extension of the continental Challenger Plateau. Depth contours are plotted from 1000-5000 m in 1000 m increments. (Norris & Cooper 2001;Sutherland et al. 2006b). Detailed bathymetric and seismic-reflection surveys demonstrate that the Alpine Fault can be traced offshore from Milford Sound to the southwest where it aligns with the southeastern edge of Resolution Ridge at the frontal Puysegur subduction thrust ( Fig. 1; Delteil et al. 1996; Barnes et al. 2005). The slip rate on the offshore Alpine Fault near Doubtful Sound is determined from offset submarine fan deposits to be 31 ± 3 mm yr -1 (Barnes 2009), which is almost all the strikeslip component of the 38 ± 1 mm yr -1 plate motion that is directed along azimuth 058 • (DeMets et al. 2010). The residual plate motion is accommodated at the surface by reverse faulting and folding of basinal sediments at 1-5 mm yr -1 west of the Alpine Fault Wood et al. 2000; Barnes et al. 2002) and brittle faulting of the upper trench slope and onshore Fiordland Sutherland et al. 2006aSutherland et al. , 2009). Present-day campaign GPS velocities indicate similarly high rates of long-term deformation on offshore structures including the Alpine Fault (Wallace et al. 2007).
To the west of the Puysegur subduction thrust, rock dredges, gravity and magnetic data show that Paleozoic continental crust of Resolution Ridge adjoins Eocene ocean crust to the southeast (Sutherland 1995;Wood et al. 1996; Barker et al. 2008). Resolution Ridge is the most southern projection of the much larger continent of western New Zealand and the Challenger Plateau, which was rifted from Australia in Late Cretaceous time and then from the Campbell Plateau during Eocene time (Sutherland 1995;Wood et al. 1996;Barker et al. 2008). It has been proposed that this Eocene rift boundary is the underlying reason for why the Alpine Fault exists: continental crust of the Australian Plate is too buoyant to subduct, so the inherited boundary with Eocene and Oligocene oceanic crust has necessarily become the locus of active deformation that connects slip on the subduction interface with the up-dip intracontinental plate boundary . As the Australian Plate moved northeast, South Tasman ocean crust was subducted at the Puysegur trench while oblique continent-continent collision took place northwest of the Eocene crustal boundary. As complexities of the Eocene rift margin were incorporated into the Neogene plate boundary, asperities were removed by displacement on the Alpine Fault and this is why the surface trace of the Alpine Fault is most complex at its southwestern end (Barnes et al. 2005). There is some debate as to whether this Eocene rift boundary has itself been reactivated to create a deep tear in the subducted Australian Plate that detaches oceanic from continental lithosphere Malservisi et al. 2003), or, as we believe is more likely, the oceanic plate is translated subhorizontally to the northeast as it is bent vertical and passes the Fiordland restraining bend Reyners et al. 2002).

1267
The Wadati-Benioff zone beneath Fiordland is entirely related to subducted Eocene-Miocene ocean crust of the South Tasman Sea Reyners et al. 2002). The architecture of the Puysegur subduction zone at depth resembles an inverted ploughshare (Christoffel & Van der Linden 1972) as the subducting plate bends around a high velocity body in the mantle (Reyners et al. 2002). Earthquake tomographic studies and hypocentres image a subducted slab that is bent to vertical at 100 km depth beneath northern Fiordland (45 • S), but is less tightly curved beneath southern Fiordland (46 • S) (Eberhart-Phillips & Reyners 2001;Reyners et al. 2002). Studies of records of large earthquakes imply a complex pattern of plate motion partitioning and internal deformation either side of the subduction thrust and Alpine Fault (Anderson et al. 1993;Doser et al. 1999;Moore et al. 2000;Reyners & Webb 2002;Robinson et al. 2003;Reyners et al. 2003;McGinty & Robinson 2007).
In this paper, we describe GPS and differential interferometric synthetic aperture radar (DInSAR) observations that define the 2009 coseismic and early post-seismic ground displacements, and we invert these observations to determine slip at the subduction interface. The 2009 earthquake occurred due east of Resolution Ridge and northeast of the Puysegur trench, and slip was mainly between Australian Plate Eocene-Oligocene ocean crust and the Pacific Plate above. The slip distribution that we derive substantially improves our understanding of how motion is accommodated at depth across the plate interface and how plate boundary motion is partitioned in the region, and we consider the implications of this event for possible future earthquakes in the region.

DATA
All continuous GPS (cGPS) stations in the southern half of the South Island recorded displacements at the time of the earthquake, with the movement exceeding 300 mm at the nearest station (PYGR) and exceeding ∼10 mm at eight stations as far away as Haast (HAAS) and Waimate (WAIM) (Figs 2 and 3). Only PYGR is in the near field of the earthquake, so the cGPS data do not provide significant constraints on the slip distribution. They do, however, provide some constraint on the overall magnitude of the earthquake.
A number of previously-occupied campaign GPS sites exist in the vicinity of the earthquake and 27 of these were reoccupied over a 5-day period 5 weeks after the earthquake (Fig. 2).
Satellite radar images have been taken over the earthquake region by the Synthetic Aperture Radar (SAR) instrument (PALSAR) on board the Japanese Space Agency's Advanced Land Observing Satellite (ALOS). The L-band (∼24 cm wavelength) radar on this satellite has been demonstrated to provide substantially more coherent images than C-band (∼5 cm wavelength) over forested and vegetated regions (e.g. Samsonov & Tiampo 2010;Zhang et al. 2010). We have processed both ascending and descending data from postearthquake scenes and combined them with pre-earthquake scenes to provide images of ground deformation during the earthquake (Fig. 4), using the DInSAR technique (Massonnet & Feigl 1998;Rosen et al. 2000). The first two post-earthquake data sets were obtained within about a week after the earthquake and the third 6-weeks later ( Table 1). The DInSAR system measures displacement along the line of sight from the ground to the satellite, so the measurement contains a mix of vertical displacement together with horizontal displacement along the azimuth from the ground to the satellite. The radar beam from PALSAR has an angle of incidence of ∼39 • , so the instrument is about 20 per cent more sensitive to vertical than to horizontal displacement.
DInSAR images (interferograms) are subject to a number of noise sources, in particular the following four: (1) Long-wavelength errors occur if the satellite orbit is not known precisely.
(2) The technique requires a correction for topography and errors arise if the digital elevation model used in the processing is not exact. We have used the 40 m digital elevation model (DEM) provided by Land Information New Zealand (LINZ). The topographic errors are proportional to the 'perpendicular baseline', or distance between the positions of the satellite when the before and after images are taken. For the images we have processed so far the baselines are between about . Time-series of daily solutions for the east and north components of the seven cGPS stations closest to the earthquake. The time-series have been regionally filtered and have had a linear trend removed so that the pre-earthquake data have zero slope. The westward and (to a lesser exent) southward displacements (except at BLUF, which moves slightly north) are clearly visible, as is post-seismic motion that continues for at least several weeks at some sites. Note that the vertical scale for the north component is double that for the east component.
400-1350 m, so any topographic error from this source should be small to moderate. (3) The radar signals are delayed as they pass through the ionosphere and through water vapour in the troposphere and these delays are indistinguishable from true ground displacement when only a pair of images is used to form an interferogram. Tropospheric errors are typically no larger than ∼10 cm and they can be correlated with topography under certain weather conditions. L-band radar is 16 times more sensitive to ionospheric electron density than C-band, so ionospheric errors can be a significant problem.
(4) The nature of the ground surface causes changes in the phase of the reflected radar beam, which is indistinguishable from a true ground displacement in a single interferogram. This is of some concern for the Dusky Sound images, as some of the 'before' and 'after' images have been taken at different times of the year, when snow depths in the more elevated regions of Fiordland are likely to be different. However, given the amount of ground deformation caused by the earthquake, all of these noise sources are likely to be small compared to the size of the coseismic signal. The images we use are summarised in Table 1. As well as the two scenes in Fig. 4 (paths 349 and 639), we use data from an ascending path (path 348) to the east of Fig. 4(a) that overlaps with the Fig. 4(a) image. This image does not contribute greatly to the solutions so it is not shown here, but all images used are reproduced in the Supplementary Information.

GPS data analysis
We analyse both the cGPS and campaign GPS data with Bernese v5.0 software. IGS final orbits and earth orientation parameters are held fixed. We use relative (rather than absolute) IGS antenna phase patterns; the GPS data used in this paper are collected almost exclusively with Trimble Zephyr Geodetic antennas so there are no problems arising from antenna mixing. We use ocean load estimates from Topex model 7.1 as calculated by Onsala Space Observatory. The tropospheric zenith delay is estimated hourly, with tropospheric tilt estimated daily, using Niell mapping functions. Ambiguities are fixed as far as possible using the Bernese software's quasi-ionosphere-free strategy and the global ionosphere model produced daily by the Centre for Orbit Determination in Europe. Daily solutions are transformed to the ITRF2000/IGb00 reference frame using a 3-parameter Helmert transformation onto a set of regional (Australia and Pacific) IGS stations. For a GPS network of the spatial extent being analysed here there is no advantage to using more recent ITRF2005/IGS05 orbits or to using absolute antenna phase patterns.

GPS displacement estimation
The displacements at the cGPS sites are estimated from regionally filtered (Wdowinski et al. 1997;Beavan 2005) daily position timeseries solutions at three post-earthquake epochs ( Table 2). The preearthquake data are the average of the solutions for July 12-14. The three post-earthquake epochs are: (1) the last 14 hours of data from July 15 (to correspond as closely as possible with the first ascending DInSAR image); (2) the average of solutions for July 22-24 (to correspond with the first descending DInSAR image) and (3) the average of solutions for August 19-23 (to correspond with the time of the GPS campaign). We take the difference between the pre-and post-earthquake observations and assign errors as the rms scatter of the before and after measurements added in quadrature. It can be seen from Fig. 3 that significant post-seismic motion occurred at some sites, most notably on the north component of PYGR. We expect that this is due to continued afterslip on the fault rupture surface spreading southward from the main rupture area, but triggered slow slip on other nearby faults could also be an explanation.
The displacements at campaign GPS sites are harder to estimate because of the several-year gap between the pre-and postearthquake observations. In the Fiordland region, we have several pre-earthquake observations at most sites that we can use to estimate the pre-earthquake motion of the site. However, many of these sites have been displaced by previous earthquakes, notably the 2003 M W 7.2 Secretary Island earthquake (Reyners et al. 2003), the 2004 M W 8.1 Macquarie earthquake (Hayes et al. 2009;Watson et al. 2010) and the 2007 M W 6.7 Fiordland earthquake. We have made corrections to the pre-earthquake station positions based on our own dislocation models of these earthquakes, which we document in Tables S1-S4 of the Supporting Information. The largest correction for the 2003 earthquake is ∼90 mm at site DF5D (Fig. 2), for 2004 it is ∼12 mm at B03W and for 2007 it is ∼20 mm at A1TH. We then estimate the 2009 coseismic displacements by fitting a straight line to the pre-earthquake data, as shown by two examples in Fig. 5. For each site we project the line and the uncertainty of the fit forward to the time of the post-earthquake GPS campaign, then subtract the measured position from the projected position to give the displacements shown in the final section of Table 2 (see also Fig. 8). The maximum observed horizontal displacement is over 800 mm and the maximum vertical is more than 200 mm subsidence, both in the vicinity of Resolution Island.    Note. a cGPS displacement estimates are made at three times following the earthquake, to correspond approximately to the times that post-earthquake DInSAR images were acquired and campaign GPS data were collected. These times are

DInSAR analysis
We analysed the SAR data with Gamma software (Wegmüller & Werner 1997) using standard DInSAR processing from ALOSspecific raw format L1.0. Adjacent frames were concatenated in one strip and processed simultaneously. A 5×10 multilook was performed to give a resolution of 23×33 m and the topographic signal was removed using the LINZ 40 m DEM. Adaptive filtering (Goldstein & Werner 1998) was performed in order to smooth the observed phase and the Minimum Cost Flow algorithm (Costantini 1998) was used for phase unwrapping. Orbits were not corrected because no significant orbital ramps were observed in the calculated interferograms; however, various linear and quadratic corrections were tested and applied during the inversion (Section 4).
Information on the images is provided in Table 1. We show examples of ascending and descending images in Fig. 4. Taken together these indicate a mixture of ground subsidence and generally westward horizontal ground displacement, which is qualitatively similar to the campaign GPS observations.

M O D E L L I N G
We invert the cGPS, campaign GPS and DInSAR data together to estimate the slip distribution and magnitude on the fault plane, which we assume from initial geodetic and seismic models (Fry et al. 2010) to be the subduction interface between the descending Australian Plate and the overlying Pacific Plate. We use a profile of the subduction interface through the Dusky Sound epicentre as estimated from the limited relocated microseismicity currently available. This profile was drawn prior to accurate knowledge of the 2009 earthquake location and therefore was not constrained to pass through the 2009 hypocentre (Fig. 6, M.E. Reyners, personal communication, 2009; M.E.R notes that the definition of the plate interface in this region will be greatly improved when the aftershocks of the 2009 earthquake have been fully processed and relocated). We form the model fault surface by extending this profile 120 km along strike to fully overlap the aftershock zone of the earthquake. We divide Figure 6. Subduction interface surface near the epicentre estimated from microseismicity relocated with a 3-D velocity model (Reyners, personal communication, 2009;Fry et al. 2010). We solve for slip between 5 and 50 km depth (red curve) and 120 km along strike, assuming the interface geometry remains the same along that length. See Fig. 8(a) for the surface projection of this area. the fault surface between 5 and 50 km depth into 5-km square cells and solve for the slip in each cell. The geographic location of the assumed fault surface can be seen by referring forward to Fig. 8. We use linear inversion, closely following the methods adopted by Jónsson et al. (2002). We use Laplacian smoothing to stabilise the solution, choosing the weighting parameter for an optimum trade-off between misfit and solution roughness. We weight the slip magnitude towards zero at the lateral and lower boundaries of the fault surface to damp artefacts that sometimes appear at the boundaries, but put no constraints on the upper boundary. This is because we wish to test if the data can constrain the upper depth of significant slip. We also tested a non-linear least-squares algorithm restricting the slip to have no normal-faulting or left-lateral components, but found that this makes only minor differences to our models.
Each interferogram contains millions of pixels, with nearby regions of the interferogram being highly correlated with each other and thus containing no independent information. We average the individual pixels into larger pixels of linear dimension ∼100 m, then apply quadtree partitioning (e.g. Jónsson et al. 2002) to further reduce the number of data points while retaining the statistically significant part of the signal. The intention of the quadtree partitioning, which we show in Supplementary Figs S1 and S2, is to obtain a set of statistically independent values that contain the important information from the original image. We assign uncertainties to the data points equal to 10/ √ A mm, where A is the area in km 2 of the region contributing to the data point. The factor 10 is chosen so that the uncertainty is 10 mm for data averaged over a 1-km square and the square root assumes the scatter about the mean in each region follows a white-noise model. Later we partially test the assigned errors by comparing the misfit of the DInSAR and GPS data sets to the model. The resulting data for inversion consist of: (1) 27 3-D displacements from campaign GPS; (2) 14 3-D displacements from cGPS (most of these are in the far-field so do not contribute much to the variable slip solution); (3) between 50 and 450 points in each of the DInSAR images, with the number depending on the smoothness of the original image.
In the modelling, we use the precise look angle from the ground to the satellite (which varies by a few degrees across the image) and we solve for an offset and a planar slope of each interferogram. The offset estimate is essential as we do not know the zero point on the interferogram, while the plane corrects for first-order orbit errors. We have also experimented with fitting a quadratic surface to correct for higher-order orbit errors. For the coseismic images this is probably not appropriate because the earthquake deformation covers much of each image and some part of the true signal may be removed by a quadratic correction. For the post-seismic images, where the much smaller deformation should tend to zero away from the earthquake source, we have obtained solutions using both planar and quadratic corrections.
Post-earthquake DInSAR data were collected on July 15 (only 2 hr after the earthquake and possibly while substantial afterslip was still occurring), July 23 and over the interval August 19-September 9 which also corresponds to the August 19-23 period of campaign GPS data collection. We refer to July 15 as epoch 1, July 23 as epoch 2 and August 21 as epoch 3. We call our models for these three epochs CP1 (Coseismic + Post-seismic at epoch 1), CP2 and CP3. We also use the cGPS and DInSAR displacements between July 15 and August 21 to determine the early post-seismic deformation. We call this model PP (Post-seismic, Planar) when we solve for a planar surface in the DInSAR data and PQ (Post-seismic, Quadratic) when we solve for a quadratic surface.
The most uniform set of post-earthquake data in terms of time of collection consists of the campaign GPS collected on August 19-23, the cGPS offsets estimated at the same epoch, the ascending images collected on August 21 and 30 and the descending image collected on September 9. These contain the coseismic deformation and the first ∼6 weeks of post-seismic deformation. Based on the cGPS time-series (Fig. 3) we do not expect large postseismic changes over the August 19-September 9 time interval of this data set; the position of PYGR changes by <5 mm, or <1.5 per cent of the signal, over this period. We first invert this data set (model CP3) to determine a slip distribution for the sum of the coseismic and first several weeks of post-seismic deformation. We then try other analyses using data from closer to the time of the earthquake to attempt to determine the evolution in postseismic slip. Table 3 summarises the various data sets we have modelled.
The equations to be solved may be written in matrix form as d = Gm where d is the data vector of observed ground displacements, m is the vector of model parameters (strike-slip and thrust displacements for each fault patch) and G contains the data kernels relating unit slip on a fault patch to surface displacement at an observation point. The data require weighting so that data points with lower uncertainty are given more weight in the solution. The weighted system of equations can be written d = G m, where d = Wd, G = WG and W is the weight matrix from Cholesky decomposition of the data covariance matrix (e.g. Jónsson et al. 2002). In our case, where we have taken each data point to be independent, the weight matrix is diagonal with elements equal to 1/σ i , where σ i is the standard error of the ith data point. Before inversion we add two blocks at the bottom of the weight matrix and data vector, as Here, D is the Laplacian of the fault slip and these terms cause the solution to become smoother as the value of the Lagrange multiplier κ 2 is increased. B is an operator that causes slip to tend towards zero at the side and lower boundaries of the model fault plane, with the effect being greater as λ 2 is increased.
Our strategy for the modelling is to give high weight to data sets collected at the epoch of interest and to either not use or give low weight to data sets collected at other times (Table 3). Ideally one would not use data sets collected at different times, which would mean not using the campaign GPS at all except at epoch 3. However, the campaign data set provides values of ground deformation that are both 3-D and absolute (though with some uncertainty attached). These are not available from the DInSAR, in the first case because each DInSAR image gives only 1-D data along the line of sight and in the second case because there is no stable reference area in the DInSAR images because the deformation covers most or all of each image. We have therefore included the campaign data at low weight at epochs 1 and 2 to provide extra stability to the inversion. We have chosen a weighting of 1.0 for the cGPS data set in all analyses. We have then adjusted the weight for the other data sets collected at that epoch (the 'primary' data sets) so that the χ 2 per degree of freedom (which we abbreviate as χ 2 from here onwards) for each of these data sets matches the χ 2 of the cGPS data set within a factor of 2 (Table 3). We then choose much lower weights for any other data set included in the inversion so that their χ 2 values are small, typically 20 per cent or smaller than the χ 2 values for the primary data sets.
To provide this additional weighting we multiply the appropriate rows of the weight matrix, W, by the weights given in Table 3.

R E S O L U T I O N, A C C U R A C Y A N D S M O O T H I N G
The data distribution for the inversion is not ideal because the data only exist onshore, whereas the slip is located partially offshore to the west. In addition, if slip has occurred well south of the epicentre, as determined by preliminary models (Fry et al. 2010), even onshore data are lacking close to the slipping region. To test the resolution we input a checkerboard slip pattern and see how well it is reproduced by the inversion. Fig. 7(a) shows an input pattern consisting of 30 × 30 km patches alternating between 4 m thrusting and zero slip. To generate the test data we forward model the 3-D displacement at each GPS station and the line-of-sight displacement at each quadtree-partitioned DInSAR point. We perturb each test data value by adding an amount randomly selected from a Gaussian distribution with a width equal to one to two times the standard error of the actual data value at that location. We then invert the test data to determine slip on the fault surface, using the same methodology as for the actual data. One example of such an inversion is shown in Fig. 7(b). From running a number of cases, we find that the method does not well resolve the checkerboard if its linear dimension is 20 km or smaller. For the 30 km case we see that the resolution is quite good in the epicentral region, with the slip amplitude determined at about the 20 per cent level (though smoothed compared to the input values). Further up dip and especially further southwest, the resolution deteriorates. We find similar results when we input a strike-slip checkerboard pattern.
We use two hyper-parameters, κ 2 and λ 2 to control the degree of smoothing of the slip distribution and the trend towards zero slip at the side and lower boundaries of the fault plane. The value of κ 2 is chosen by a commonly used method of plotting data misfit against solution roughness as κ 2 varies between 3000 and 10 ( Fig.  S3 of Supporting Information) and choosing a value where the misfit is falling slowly but the slip distribution is reasonably smooth (e.g. Jónsson et al. 2002). We use the same value, κ 2 = 100, for all the coseismic data sets we invert. For the post-seismic data sets we use higher values (Table 3) to give smoother solutions. In Figs S4-S9 of the Supporting Information we investigate the sensitivity of the model solution to (1) different values of the hyperparameters (using smaller and larger values of κ 2 and setting λ 2 to 0), (2) increasing the area of the model fault surface and (3) changing the relative weighting of the GPS and DInSAR data. These tests give us confidence in our model results. Figs 8,9 and 11(c) show the inversion results and the fit to the data for model CP3, where we have weighted the individual data sets so that each has approximately the same reduced χ 2 ( Table 3). The majority of slip occurs in the region where the model resolution is good (Figs 7 and 11c). In particular, several features are well defined: the fairly sharp cut-off of slip at the NNE end of the rupture; the ∼5 m maximum slip centred at ∼17 km depth about 30 km SW of the epicentre; the maximum depth of significant slip at about 35-40 km depth in the epicentral region and the reduction of slip at shallow depth towards the upper end of the fault surface. Fig. 11 shows the slip distribution on the subduction interface for the various models. Most slip has already occurred by the time of model CP1 (Fig. 11a, Table 3). Given the resolution test of Fig. 7 it is hard to be sure whether the differences between CP1, CP2 and CP3 (Figs 11a-c) are significant. However, the cGPS data clearly show that post-seismic deformation is occurring and the interferograms between the post-seismic July and August/September images do suggest a signal associated with the quake, though orbital errors may also contribute to the signal (Fig. 10).

R E S U LT S
Models PP and PQ (Figs 10, 11d and e) attempt to isolate the afterslip pattern that caused the July 15-August 21 post-seismic deformation. If all the data were noise-free and the inversions perfect, then models PP and PQ (Figs 11d and e) should give identical slip distributions, which in turn should be the same as the difference between the CP3 and CP1 slip distributions, which we show in Fig. 11(f). There are some similarities: there is post-seismic slip below and/or to the right (NNE) of the epicentre in all cases and there is post-seismic slip to the left (SSW) of the main slip patch. However, in models PP and PQ the slip patch to the SSW is 5-10 km deeper than in the difference between CP3 and CP1. We note that modelling the difference between two observed data sets is a more stable procedure than taking the difference between two independently-calculated models, so we prefer the post-seismic models of Figs 11(d) and (e) over that of Fig. 11(f).
The inferred moments and maximum slip estimates increase approximately monotonically from model CP1 to CP3, which is consistent with true post-seismic deformation being recorded in the data. Also the moments from models PP and PQ are only about 25 percent greater than the difference between models CP3 and CP1 (Table 3) which gives further confidence that the DInSAR is detecting post-seismic motion even though the details are not very well resolved.
In the Supplementary Information (Figs S10-S13) we show similar DInSAR images to those in Figs 9-10, but for models CP1, CP2 and PP and for the track 348 image used in model CP3.

Model verification
The observed and modelled subsidence in the vicinity of Resolution Island and Dusky Sound are both in the 200-250 mm range. Initial field surveys in the region for coastal uplift and subsidence detected inundation of salt marsh plants indicating the possibility of 100-150 mm subsidence, with 300 mm the maximum possible subsidence compatible with the observations ). These observations are in good agreement with the geodetic data and models.
Tsunami waveforms have been predicted by Prasetya et al. (2010) using a coseismic geodetic model very similar to the ones described here. The model matches the observed waveforms at DART buoy 55015 in the Tasman Sea very well without any additional tuning, which also supports the accuracy of our coseismic displacement model.

Model fit to data
None of our model fits achieve a χ 2 of 1. The best we achieve is χ 2 ≈ 2 for model CP1 and χ 2 ≈ 3-4 for models CP3, PP and PQ. This suggests that our data uncertainties are underestimated by factors of 40-100 per cent and/or that there are systematic inconsistencies between the data types. This level of fit is, however, similar to that observed in many other studies of coseismic slip using geodetic data. The fit for model CP2 is worse, with χ 2 ≈ 8-11 when we balance the fit to the cGPS and DInSAR data sets. The reason for this may be that the DInSAR descending data for 2009 July 23 contains more noise signals than the similar data for August 30. Fig. S5 in the Supporting Information shows a solution using the same data and hyper-parameters as model CP3 but with the DInSAR data down-weighted by a factor of ∼30 compared to CP3. A comparison of Fig. 11(c) (which is also redrawn as Fig. S4) and Fig. S5 shows some small variations in slip direction and magnitude, with the location of the major slip patch quite similar in the two cases. This indicates that the GPS and DInSAR data contain similar but not identical information about the ground deformation. We prefer the solution that combines the GPS and DInSAR data sets with comparable χ 2 values, as we have no way to know which is the more reliable data set.
This comparison also indicates that for model CP3 the DInSAR data add only a small amount of information that is not already available in the combined continuous and campaign GPS data sets. This might at first sight be seen as an argument against the need for DInSAR data. However, it is the DInSAR data, with assistance from the sparse continuous GPS, that defines models CP1 and CP2 and thus gives us what information we have about the early postseismic processes. Also, in this particular example the majority of the deformation signal is offshore and thus unavailable to either GPS or DInSAR; the relative advantages of DInSAR in measuring the near-field ground deformation at high spatial resolution are thus less apparent in this case than for earthquakes that are primarily on land. Finally, this is an example where a reasonably dense campaign GPS network is available; in regions where this is not the case, DInSAR is still able to provide a ground deformation measurement.

Coseismic slip and implications for regional tectonics
The 2009 earthquake provides the first unequivocal evidence for slip on the Puysegur subduction interface and can thus be used to test previous hypotheses for how slip partitioning is occurring at the plate boundary adjacent to Fiordland Sutherland et al. 2000;Reyners et al. 2002).  . Observed (top), modelled (middle) and residual (lower) DInSAR images for model CP3, for the August 30 ascending path (track 349, left) and the September 7 descending path (track 639, right). The images show the observed interferogram after reconstruction from the several hundred points in the quadtree decomposition, which is why the top images look slightly different to those in Fig. 4. Note that the colour range is 600-700 mm on the observed and model images, ∼200 mm on the ascending residual and only ∼60 mm on the descending residual. The poorer fit to the ascending image may be because the baseline length is significantly longer than in the descending image, so that coherence is lower and phase unwrapping followed by interpolation is less accurate.   (track 639, right). It is likely that there is a significant orbital ramp in the observed images, which has been estimated along with quadratic terms in the inversion. The presence of these relatively large orbital terms makes it difficult to be confident of the inversion for post-seismic slip (e.g. Figs 11d-f below).  The region that we model to have a slip of more than 3 m occurred at 10-30 km depth on an interface that dips 25-35 • (Fig.  6). Our best-fitting model of slip direction has a rake of 141±6 • within this zone of maximum slip (Fig. 8). The expected rake of the plate motion on this interface, accounting for rotation to this plane (bending of the slab), is 146 ± 1 • for the MORVEL model (DeMets et al. 2010) or 142 ± 2 • from regional GPS data (Beavan et al. 2002). Therefore, the slip direction during the 2009 earthquake was, within uncertainty, the same as that expected if all plate motion were occurring on the subduction interface. If all plate motion were occurring on the interface, the maximum model slip of ∼5.4 m represents ∼140 years of accumulated plate motion.
The direction of slip in the 2009 earthquake is consistent with no slip partitioning on the plate boundary at 10-30 km depth beneath Dusky Sound, but this is in very obvious contrast to the situation 100 km farther north near Doubtful Sound. Seabed mapping of fault traces Barnes et al. 2002Barnes et al. , 2005 and earthquake focal mechanisms (Anderson et al. 1993;Moore et al. 2000;Reyners & Webb 2002;Reyners et al. 2003;Robinson et al. 2003;McGinty & Robinson 2007) reveal that shallow crustal deformation near Doubtful Sound is partitioned between strike-slip motion on a steeply dipping Alpine Fault and dip-slip thrust or reverse faulting on the subduction interface. The distribution of slip that we calculate during the 2009 earthquake is illustrated in Fig. 12 and provides an explanation for this significant spatial change in deformation style. Slip during the 2009 earthquake was on the interface between Eocene-Miocene ocean crust and the overlying Pacific Plate and was (within uncertainty) entirely down-dip of the Eocene rift boundary in the Australian Plate. This inherited Eocene crustal boundary on the surface of the subducted plate plunges to lower crustal depths (ca. 30 km) near Doubtful Sound, so crustal faulting in that region is entirely within an oblique continental collision zone, even though slip on the subduction interface is presumably not partitioned at depths below where the Alpine Fault meets the subduction interface.
The geometry of deformation that we propose in Fig. 12 requires northwestwards advection of continental crust and the Alpine Fault, which is the boundary between Australian Plate and Pacific Plate continental crust, relative to the Australian Plate. This is required because shortening of Australian Plate crust is known to occur at the plate boundary. Barnes et al. (2002) estimate ca. 7 km of shortening within the Fiordland basin near Dusky Sound since 3 Ma. The Eocene crustal boundary within the Australian Plate beneath the subduction interface must, therefore, be progressively offset downdip relative to the surface trace of the Alpine Fault, and this may lead to a smearing of the Alpine Fault at depth to moderate southeast dips. This is in contrast to the subvertical plane mapped at the surface (Barnes et al. 2005), but is consistent with the distribution of dip-slip focal mechanisms and specifically the 2003 Secretary Island event, which had a dip-slip down-dip rupture extent that was significantly southeast of the Alpine Fault surface trace (Reyners et al. 2003;McGinty & Robinson 2007).
Our preferred model has the up-dip extent of the 2009 rupture approximately beneath the surface trace of the Alpine Fault, but we cannot rule out the possibility that rupture continued up-dip into Australian Plate crust on fault splays that meet the surface within the Fiordland Basin (Barnes et al. 2002). We do not require strike-slip motion on the Alpine Fault to explain our GPS observations, but we ran forward models of right-lateral slip on the Alpine Fault to see if the predicted signal would be observable at our GPS stations. To generate a displacement signal larger than the observed differences between observations and model (Fig. 8) would require several metres of slip over tens of kilometres fault length, so we cannot rule out small amounts of slip over a limited region (e.g. 1 m over a 10 km length) of the Alpine Fault. Preliminary aftershock studies do not reveal strike-slip mechanisms or greater hypocentre density near to the Alpine Fault. The Alpine Fault surface trace has not been reobserved since the 2009 rupture.
Was any secondary surface faulting triggered in Fiordland during the 2009 earthquake? We examined by eye the differences between observed interferograms (before down-sampling) and the models, to see if any 'fault-like' signals are evident in the residuals. Our characterization of the resulting images is that we see noise signals with amplitudes up to about ±10 cm at various spatial scales (some of which are visible in the residual plots of Figs 9-10), but we found no signals that are clearly indicative of on-land crustal faulting.

Post-seismic slip
With only one cGPS station in the near field it is not possible to place strong constraints on the DInSAR images to solve for a reliable history of post-seismic deformation from July 15-July 23-August 21. However, our models PP and PQ and the differences between models CP1 and CP3 do provide a somewhat consistent view of the afterslip that caused the observed post-seismic deformation. We assume that the majority of post-seismic deformation in the first few weeks is due to afterslip as has been widely observed elsewhere, but expect that ongoing post-seismic deformation may be due to longer-term viscoelastic or other effects (e.g. Freed 2005). Though the details are not well resolved, the inferred afterslip is located around the edges of the main coseismic slip patch, particularly near its southwest end and below and northeast of the epicentre. The inferred afterslip is also located in the same regions as the majority of aftershocks, suggesting that both are a response to stress changes caused by the earthquake. Barnes (2009) has established that the offshore Alpine Fault has a long-term average slip rate of ∼30 mm yr -1 at least as far south as Doubtful Sound. Assuming that the Alpine Fault continues at ∼30 mm yr -1 south of Doubtful Sound (Barnes, personal communication 2010, suspects it does because of its strong geomorphic expression, though no offsets have been dated south of Doubtful Sound) there is potential for a future earthquake on this part of the fault. Fry et al. (2010) used a simplified preliminary model of the slip distribution to show that the Coulomb failure stress increases on parts of the Alpine Fault. We have repeated that calculation using a close approximation to the CP3 slip distribution and the results are presented in Fig. 13. The Coulomb stress for right-lateral failure on the Alpine Fault increases along a ∼100 km length of the fault between Resolution Ridge and south of Doubtful Sound by an average of 0.55 MPa in the 0-15 km depth range and 0.68 MPa in the 0-10 km range that is likely to be relevant towards the southern end of the fault. This ∼100 km section of the fault has therefore been brought closer to failure; however, it is not known when an earthquake last occurred on this fault section.

C O N C L U S I O N S
We have shown that L-band DInSAR is a valuable technique for measuring coseismic ground displacement, even in a region like Fiordland that has rugged and steep topography with substantial areas of dense vegetation. Together with GPS data, differential InSAR provides a powerful method to invert for coseismic slip distributions even when ground deformation data are available on only one side of the earthquake.
Our results demonstrate oblique subduction in action. In the region of the 2009 earthquake nearly all the plate motion below 10 km depth is taken up by oblique slip on the plate interface. Only at depths shallower than 10 km is slip partitioned between strike-slip motion on the right-lateral Alpine Fault and a few mm yr -1 shortening across the Fiordland Basin to its west. This is in sharp contrast to the situation 100 km further north where slip is partitioned between strike-slip on the Alpine Fault and nearly pure thrusting on the subduction interface to depths of at least ∼25 km. We explain the difference as a result of evolution of the Alpine Fault as Challenger Plateau continental crust on the Australian Plate moves north and resists subduction beneath the continental rocks of New Zealand. If the 2009 earthquake is a typical example of how plate motion is accommodated in this region, we would expect similar events on average about every 150 yr.
We observe post-seismic surface displacements in the 5 weeks following the earthquake. Assuming that these are due to afterslip on the main fault surface we find (though with low resolution) that the afterslip occurs in regions surrounding the main slip patch, which also coincides with the majority of aftershock activity. This pattern is similar to that observed in a number of other cases worldwide.
The coseimic deformation has caused an increase in Coulomb stress promoting right-lateral strike-slip failure on a ∼100-km long section of the offshore Alpine Fault. The increase averages ∼0.6 MPa, but it is not known when this fault last ruptured nor whether rupture typically initiates on this section of the fault.

A C K N O W L E D G M E N T S
We thank GeoNet (www.geonet.org.nz), the New Zealand Earthquake Commission and Land Information New Zealand (LINZ) for the operation and financial support of the cGPS network, LINZ for participation in and partial support of the post-earthquake GPS campaign, Chris Pearson for his role in establishing the Fiordland GPS nework, Phil Barnes for discussions and Laura Wallace and Martin Reyners for their reviews of the manuscript. Also Jack Loveless, an anonymous reviewer and editor Duncan Agnew for their useful comments. The ALOS PALSAR data has been used in this work with the permission of JAXA and METI and the Commonwealth of Australia (Geoscience Australia) ('the Commonwealth'). JAXA, METI and the Commonwealth have not evaluated the data as altered and incorporated within this work and therefore give no warranty regarding its accuracy, completeness, currency or suitability for any particular purpose. Fig. 13 was created using Russell Robinson's GNStress 2.17 software, available at ftp://ftp.gns.cri.nz/pub/robinson/GNStress2. The inversion calculations were done with Igor Pro (www.wavemetrics.com), which was also used to prepare many of the figures. Other figures were prepared with GMT (gmt.soest.hawaii.edu). The majority of this work was completed with funding from the New Zealand Foundation for Research, Science and Technology.

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 article: Tables S1-S4. These give the parameters of the dislocation models and the resulting displacement values that we used to correct the campaign GPS time-series for coseismic offsets due to the 2003 Secretary Island (Fiordland), 2004 Macquarie and 2007 Fiordland (George Sound) earthquakes. These models provide a reasonable to good fit to observed data, but they are not necessarily the best dislocation model approximations obtainable for these earthquakes. Tables S5-S7. These give the parameters of our variableslip dislocation solutions CP1, CP2 and CP3 for the 2009 Dusky Sound earthquake. These are provided in separate tabdelimited text files: Beavan_DuskySound_TableS5_CP1.txt, Bea-van_DuskySound_TableS6_CP2.txt and Beavan_DuskySound_ TableS7_CP3.txt.
Figures S1-S2. These show the quadtree partitioning of the DIn-SAR images used in our analysis. Figure S3. This shows the trade-off curve between solution misfit and solution roughness that is used to select the value of the κ 2 hyper-parameter. Figures S4-S9. These show variations in the CP3 solution as the κ 2 and λ 2 hyper-parameters are varied, the dimensions of the fault surface are increased and the DInSAR data are down-weighted compared to the GPS data.
Figures S10-S13. These contain plots of the quadtree-sampled DInSAR images, model fits and residuals used in our analysis but not included in the main paper.
Please note: Wiley-Blackwell are 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 article.