Systematic assessment of atmospheric uncertainties for InSAR data at volcanic arcs using large-scale atmospheric models: Application to the Cascade volcanoes, Remote Sensing

Satellite Radar Interferometry (InSAR) is suited to monitoring ground deformation on the scale of volcanic arcs, providing insight into the eruptive cycle over both long and short time periods. However, these measurements are often contaminated with atmospheric artefacts caused by changes in the refractivity of the atmosphere. Here, we test the use of two large-scale atmospheric models, ERA-Interim (ERA-I) and North American Regional Reanalysis (NARR), to correct atmospheric uncertainties in InSAR data from the Cascades Volcanic Arc, United States. At Lassen Volcanic Center, we ﬁ nd that NARR reduces interferogram standard deviation in 79% of cases by an average of 22%. Using NARR, we develop a strategy to produce a priori estimates of atmospheric uncertainties on an arc-wide basis. We show that in the Cascades, the RMS variation in range change is dependent upon volcano topography and increases by 0.7 cm per kilometre of relief. We use this to estimate detection thresholds for long-term monitoringofsmallmagnitude(1cm/yr)deformationsignals,andshort-termmonitoringofgrounddeformation associated with pre-eruptive unrest. Thisnew approach ofassessing atmospheric uncertainties a priori iswidely applicable to other volcanic arcs, and provides realistic estimatesof atmospheric uncertainties suitable for use in near-real-time analysis of InSAR data during periods of volcanic unrest.

The potential for large-scale InSAR studies has increased with the launch of new satellites, e.g., Sentinel-1A (European Space Agency) and ALOS-2 (Japanese Aerospace Exploration Agency), providing imagery for InSAR globally every 12 and 14 days respectively (Hanssen & Rocca, 2009;Suzuki, Osawa, Hatooka, & Watanabe, 2009). However, despite advances in data acquisition, sources of noise continue to limit InSAR measurements of small magnitude ground deformation. The most significant of these is spatiotemporal variability in atmospheric refractivity between satellite acquisitions, resulting in atmospheric artefacts in interferograms that may mask or lead to false interpretations of ground deformation (e.g., Beauducel, Briole, & Froger, 2000;Poland & Lu, 2008). Improving the accuracy of satellitebased measurements of volcano deformation is essential for regional and global volcano monitoring strategies, as ground deformation is shown to have strong evidential links to eruption . When monitoring volcanic arcs using InSAR, there is therefore a need for a priori estimates of atmospheric uncertainties at a large number of volcanic edifices. Given their spatial coverage, large-scale atmospheric models are well suited to this task, and can be used to estimate detection thresholds on arc-wide scales. In this study we test this approach using two large-scale atmospheric models and InSAR data from the Cascades, which has atmospheric conditions broadly representative of a large number of other volcanic arcs (Foster et al., 2013).

The Cascades volcanic arc
The Cascades volcanic arc is located in the western US, spanning over 1000 km through Washington, Oregon and California. The arc consists of a chain of N-S trending stratovolcanoes, lava domes, cinder cones, and shield volcanoes, with 13 principle edifices ranging in elevation between 2400 m (Medicine Lake Volcano, CA) and over 4000 m (Mount Rainier, WA) (Fig. 1). Amongst the volcanoes lie numerous population centres, including the cities of Seattle and Portland, and major infrastructure such as Interstate 5. In addition to seismicity, geodetic measurements are an essential component of monitoring the arc (e.g., Ewert, Guffanti, & Murray, 2005), and InSAR data are ideal for cases where access is difficult (e.g. Mount Baker, WA), or where little ground-based equipment is deployed (e.g. Mount Adams, WA). A significant archive of InSAR data exists for the Cascades, including ERS-1/2 and ENVISAT data accessible via the WInSAR consortium, and ALOS data available from the Alaska Satellite Facility. Measurements have been made at Three Sisters, OR (Dzurisin, Lisowski, & Wicks, 2009;Dzurisin, Lisowski, Wicks, Poland, & Endo, 2006;Riddick & Schmidt, 2011;Wicks et al., 2002), Medicine Lake Volcano, CA (Parker, Biggs, & Lu, 2014;Poland et al., 2006), and Lassen Volcanic Center, CA (Poland, Bawden, Lisowski, & Dzurisin, 2004). However, overall, the application of InSAR data in the Cascades has been limited. Whilst this is partly due to poor coherence and a lack of significant deformation, measurements are also limited by extensive atmospheric noise, and published InSAR studies are of variable success. For example, interferograms from different satellite tracks covering the 2004-2008 eruption of Mount St. Helens show conflicting evidence of edifice inflation and deflation (Poland & Lu, 2008).

Causes and effects
Atmospheric phase delays can be split into a hydrostatic component, which is a function of pressure, and a "wet" component, dependent upon the water vapour content of the atmosphere (e.g., Bevis et al., 1992;Hanssen, 2001). Although the magnitude of the hydrostatic delay is several times larger than that of the wet delay, the atmospheric water vapour content is far more variable between satellite acquisitions and is therefore the dominant source of noise for differential SAR measurements (Zebker, Rosen, & Hensley, 1997). However, it is necessary to account for both wet and hydrostatic components to fully describe atmospheric phase delays in regions of significant topography (e.g., Doin, Lasserre, Peltzer, Cavalie, & Doubre, 2009;Elliott, Biggs, Parsons, & Wright, 2008;Jolivet et al., 2014). Atmospheric phase delays are commonly a few cm in magnitude (Hanssen, 2001), with a 20% change in relative humidity between acquisitions resulting in a 10 cm path delay (Zebker et al., 1997). The deformation signal measured by InSAR can therefore be modified or even reversed by atmospheric noise (e.g., Heleno et al., 2010).
Atmospheric artefacts observed in interferograms are divided into two types: stratified and turbulent (Hanssen, 2001). Vertical stratification of the atmosphere results in variable phase delays over low topography (large amounts of atmospheric water vapour) and small magnitude delays over higher topography (small amounts of atmospheric water vapour) (e.g., Ebmeier et al., 2013a). The effects are therefore most prevalent in regions of significant topographic relief including volcanic arcs (e.g., Webley, Wadge, & James, 2004). The resulting artefacts are correlated with elevation and appear as concentric fringes centred on the volcanic edifice as was first observed at Mount Etna, Italy (Beauducel et al., 2000;Delacourt, Briole, & Achache, 1998;Massonnet, Briole, & Arnaud, 1995). This is particularly misleading as volcano deformation also typically correlates with topography, and alternating periods of uplift and subsidence interpreted to reflect changes in magma storage have, in some cases, been shown to arise due to topographically correlated atmospheric artefacts. (e.g., Llaima volcano, Chile: Bathke, Shirzaei, & Walter, 2011;Remy et al., 2015). Unlike stratification, turbulent mixing of the atmosphere is not directly correlated with topography. Instead, artefacts typically exhibit spatial correlation over length scales of~10 km (e.g., Jónsson, Zebker, Segall, & Amelung, 2002;Lohman & Simons, 2005), although steep stratovolcanoes are commonly associated with turbulence on 1 km scales (Webley et al., 2004).
Volcanic arcs act as significant topographic barriers and we may therefore expect variations in atmospheric phase delays due to orographic effects (e.g., Doin et al., 2009;Puysségur, Michel, & Avouac, 2007;Wadge et al., 2002Wadge et al., , 2010. As the prevailing wind blows air towards the arc it may be forced to rise. As the air cools, clouds and precipitation form on the windward slope, causing lower levels of water vapour and precipitation on the leeward slope (Price, Byers, Friend, Kohler, & Price, 2013). At the Cascades volcanic arc, the arc axis is orientated perpendicular to the direction of the prevailing wind and within~150 km of the Pacific Ocean ( Fig. 1). As such, atmospheric conditions are strongly influenced by the interaction between the moist marine boundary, the drier air of the interior, and the complex orography (Foster et al., 2013). The result is a significant contrast between the environment of the windward and leeward sides of the arc. This is observable in maps of precipitation and the Köppen climate classification (Peel, Finlayson, & McMahon, 2007), one of the most widely used climate classification systems (Fig. 1). West of the Cascades, the climate is marine-dominated, with moderate temperatures and persistent winter precipitation (Mass, 2008). To the east, the climate is continental, with hot summers, cold winters, and low levels of precipitation (Price et al., 2013). Contrasts are also apparent between the north and south of the arc, as volcanoes at lower latitudes receive higher levels of solar radiation than those in the north.

Impact upon interferogram coherence
At volcanoes, snow cover and vegetation cause scatterers on the ground to vary rapidly over time, significantly decreasing the coherence of InSAR data (e.g., Lu & Freymueller, 1998;Lu, Power, McConnell, Wicks, & Dzurisin, 2002). The extents of both snow and vegetation cover relate to atmospheric conditions. Volcanoes that experience lower temperatures, due to latitude or elevation, will have greater snow cover resulting in poorer coherence. For example in the Cascades, snow cover has inhibited the use of InSAR at Mount Baker, WA, in the north of the arc. Volcanoes that are subject to orographic effects will have more vegetation on the windward flanks where precipitation is greater. For example, west to east across the Cascades arc, vegetation changes from dense forests to shrubs and grasses, significantly reducing the coherence on the western flank of Three Sisters, OR (Dzurisin et al., , 2009Riddick & Schmidt, 2011;Riddick, Schmidt, & Deligne, 2012;Wicks et al., 2002).
Volcano size also influences atmospheric conditions and the development of vegetation (the Massenerhebung effect), as the larger the mountain "mass", the greater the heat retention and the slower the rate at which temperature decreases with elevation (e.g., Bell, 2012;Grubb, 1971). At Three Sisters Volcano, OR, where three peaks join to form a relatively large land mass above 1800 m, alpine vegetation is well developed (Price et al., 2013) and we may therefore expect poor coherence. This contrasts to the less massive Mount Hood, OR, where the timberline is 150-300 m lower, alpine vegetation is considerably more impoverished (Price et al., 2013), and vegetation is therefore likely to have a lesser impact upon coherence. In combination, atmospheric artefacts and incoherence make interferogram phase-unwrapping extremely challenging on the steep slopes of stratovolcanoes (Pinel et al., 2011).

Atmospheric corrections for volcano InSAR studies
Multiple interferograms can be used to reduce the effects of atmospheric noise and improve the signal-to-noise ratio of InSAR imagery. If the atmospheric component is random in time, stacking N independent interferograms will reduce the standard deviation of atmospheric errors by a factor of ffiffiffiffi N p (e.g., Biggs, Wright, Lu, & Parsons, 2007;Emardson, Simons, & Webb, 2003). However, this is less useful for removing stratified delays, which do not vary randomly in space or time (Doin et al., 2009), or in the case of discrete deformation events that may be covered by a single or small number of interferograms. In cases of transient deformation, a time-series approach (e.g., SBAS; Berardino, Fornaro, Lanari, & Sansosti, 2002) may be more useful, but there is a risk of aliasing seasonal atmospheric effects into time-series analysis (e.g., Doin et al., 2009). As time-series measurements are made relative to a far field region of the interferogram, local differences in atmospheric delay may also dominate the time-series (Ebmeier et al., 2013a). Beyond simple time-averaging approaches, significant progress has been made towards correcting atmospheric noise in InSAR data, including a number of studies that tackle the problem in volcanic settings. These methods are broadly categorised in two groups: empirical and predictive (e.g., Jolivet et al., 2014).

Empirical atmospheric corrections
The most common empirical approach to correcting atmospheric errors is quantifying the correlation between elevation and delay to identify phase delays due to atmospheric stratification (Delacourt et al., 1998). This may be done by plotting the phase of interferogram pixels against corresponding elevations from a digital elevation model. From this, a correlation coefficient and elevation-delay relationship can be calculated. The elevation-delay relationship may be approximated as linear (e.g., Cavalié, Doin, Lasserre, & Briole, 2007;Wicks et al., 2002), although for steep-sided volcanoes a non-linear model may be more appropriate (e.g., Remy, Bonvalot, Briole, & Murakami, 2003). In the presence of lateral heterogeneities in stratification, for example due to orographic effects, a single empirical relationship may not be applicable across a whole scene (e.g., Puysségur et al., 2007). More advanced empirical methods therefore involve prediction of the elevation-delay relationship over different spatial scales using either Gaussian filters of varying width (Lin, Simons, Hetland, Muse, & DiCaprio, 2010), or wavelet transforms to identify correlations over varying spatial frequencies (Shirzaei & Bürgmann, 2012). However, as volcano deformation typically correlates with edifice topography, empirical corrections are particularly likely to introduce artefacts or remove elements of the deformation signal (Ebmeier et al., 2013a).

Predictive atmospheric corrections
Predictive methods utilise inputs from external sources to estimate atmospheric parameters and produce maps of the expected atmospheric delay. This includes the use of large-scale atmospheric models (e.g., Popocatepetl and Colima Volcano: Pinel et al., 2011), datasets such as GPS (e.g., Mt. Etna: Li, Ding, Huang, Wadge, & Zheng, 2006), local weather data (e.g., Mt. Etna: Delacourt et al., 1998), plus estimates of water vapour from multi-spectral imagery (e.g., Fogo Volcano and Mount Cameroon: Heleno et al., 2010). In addition to large-scale atmospheric models, mesoscale models have also been used in volcanic settings and are often implemented in a nested fashion within the coarser grid of a large-scale model (Eff-Darwich et al., 2012;Foster et al., 2006).
Atmospheric corrections using predictive methods have been of varied success. At Mount Etna, Li, Ding, et al. (2006) and Webley et al. (2004) report a~30% reduction in interferogram phase standard deviation following corrections using GPS data and forward models of the atmosphere respectively. However predictive methods are not always successful: in their study at Mount St. Helens, Foster et al. (2013) assimilated ground-based and remote sensing data with the predictive, mesoscale MM5 atmospheric circulation model, but found this provided no mean benefit for interferogram analysis.
The applicability of predictive methods is dependent upon the availability of external data. Multi-spectral imagery, such as that provided by MERIS and MODIS, relies upon cloud free, day-time conditions, and ground-based measurements/GPS are typically sparse, except for intensely monitored volcanoes. Many external datasets also suffer from timing issues, as data are not acquired at the same time as InSAR passes, and atmospheric conditions vary rapidly over timescales of b1 h (Hanssen, 2001). In this respect, atmospheric models are advantageous: large-scale atmospheric models provide multiple simulations each day (e.g., Pinel et al., 2011), and nested models can be run to any time step (e.g., Foster et al., 2006). However, these results are not necessarily accurate as atmospheric models are notoriously sensitive to initialisation conditions (e.g., Wadge et al., 2002) and the time at which atmospheric phenomena are modelled to have occurred may be offset from the actual time (Zhu et al., 2007). The large grid spacing of regional and global-scale models may also be insufficient to simulate smaller scale atmospheric heterogeneities such as those associated with turbulent mixing (Jolivet et al., 2014;Webley et al., 2004). So far, these predictive methods have been applied to individual volcanoes on a case-by-case basis (see references mentioned throughout this section). This study therefore presents the first application of large-scale atmospheric models to assess atmospheric uncertainties over an entire volcanic arc.

InSAR
We test the application of large-scale weather models in the Cascades using examples of InSAR data that show artefacts due to both atmospheric stratification and turbulent mixing. By assessing past deformation studies in the Cascades (discussed in Section 1.1) and processing test scenes for volcanoes throughout the arc, we select two case study volcanoes, Lassen Volcanic Center and Medicine Lake Volcano, which demonstrate atmospheric noise of varying characteristics. We choose these case study volcanoes as they are located in regions with the driest, warmest Köppen classification ( Fig. 1), and are therefore most likely to have good coherence due to minimal snow cover and vegetation. They also represent two different types of volcano topography as defined by the Smithsonian Global Volcanism Programme: Medicine Lake Volcano is classed as a broad shield volcano with summit caldera, whereas Lassen Volcanic Center is a stratovolcano. For each volcano we use a dataset summarised in Table 1, with a full list of interferograms supplied in the Supplementary Material. Deformation rates at these volcanoes during the observation period have been measured to be~1 cm/yr. (Dzurisin, Poland, & Bürgmann, 2002;Poland et al., 2006). Rather than masking the volcanic edifice (e.g., Doin et al., 2009), we use interferograms with duration b 1 year (maximum deformation signal 1 cm), which allows us to minimise any bias caused by topographically-correlated deformation signals whilst preserving the region of interest in the interferograms.
Interferograms from Lassen Volcanic Center exhibit topographically correlated phase signals and therefore present examples of noise related to atmospheric stratification. The dataset covering this volcano consists of 44 C-band SAR images from ENVISAT ascending track 435 ( Fig. 1), which are used to produce 38 interferograms spanning May 2004-August 2010. The second dataset, covering Medicine Lake Volcano, has been used in past InSAR studies and contains short wavelength phase heterogeneities that are not correlated with topographic features (Parker et al., 2014). These interferograms present examples of turbulent atmospheric phase delays that are known to be problematic for InSAR studies. The dataset covering Medicine Lake Volcano consists of 32 acquisitions from ENVISAT descending track 342, which are used to produce 54 interferograms spanning November 2005-October 2010.
All interferograms used were produced using the JPL/Caltech ROI PAC software (Rosen, Hensley, Peltzer, & Simons, 2004), and a 90 m SRTM digital elevation model (Farr et al., 2007). Interferograms were filtered using a power spectrum filter (Goldstein & Werner, 1998), unwrapped using a branch cut algorithm (Goldstein, Zebker, & Werner, 1988), downsampled to a final resolution of~500 m (e.g., Goldstein et al., 1988;Jónsson et al., 2002), and then converted from phase to millimetre range change with positive displacements corresponding to movement towards the satellite.

Atmospheric models
We use two atmospheric models, ERA-Interim (ERA-I) and the North American Regional Reanalysis (NARR). Both of these atmospheric models are based upon data reanalysis such that past observations of the atmosphere are reassessed using a numerical weather prediction forecast and analysis system (Berrisford et al., 2011). Both ERA-I and NARR are freely available online along with documentation describing the format of the data.
ERA-I is a global atmospheric model from the European Center for Medium-Range Weather Forecasts (ECMWF) (Dee et al., 2011;Uppala, Kobayashi, Berrisford, & Simmons, 2008;Simmons, Uppala, & Kobayashi, 2007). The model is produced using a sequential data assimilation scheme, the core component of which is a 4-dimensional variational analysis, which advances in 12-hourly analysis cycles. Within each cycle, observations are combined with prior information from a forecast model to estimate the evolving state of the global atmosphere (Dee et al., 2011). The number of observations assimilated each day is of the order of 10 7 , with the majority of data originating from satellites. The analyses are then used to initialise a short-range model forecast, which provides the prior state estimates needed for the next analysis cycle (Dee et al., 2011). The outputs of ERA-I used in this study are estimates of temperature, relative humidity, and geopotential height, defined at 37 pressure levels (1000-1 hPa), and a spatial resolution of 75 km (Fig. 1). NARR is from the National Center for Environmental Prediction/ National Center for Atmospheric Research (NCEP/NCAR) and provides coverage for North and Central America (Mesinger et al., 2006). The NARR assimilation system operates over 3 -hour cycles, with the preceding cycle serving as the first estimate for the next cycle (Rogers et al., 2001). Data sources include rawinsondes, dropsondes, aircraft, satellite and surface based measurements such as observed Table 1 InSAR datasets used in this study and shown in Fig. 1 precipitation (Mesinger et al., 2006). Like ERA-I, NARR defines a set of meteorological parameters, including temperature, geopotential height, and specific humidity, at 29 pressure levels (1000-100 hPa), and a spatial resolution of 32 km (Fig. 1). For both models, pressure levels are more densely spaced at lower elevations (higher pressure) and widely spaced at higher elevations (lower pressure). We select the cycle that is closest to the time of SAR acquisition, which for both models is within 30 min of ENVISAT ascending pass (~0600 UTC) and ENVISAT descending pass (~1800 UTC). To produce maps of atmospheric delay, we convert relative and specific humidities to partial pressure of water vapour, and geopotential to altitude. We then interpolate the atmospheric parameters onto altitude profiles at each model node using a spline interpolation (Jolivet, Grandin, Lasserre, M.-P., D., & Peltzer, 2011;Jolivet et al., 2014;Walters, Elliott, Li, & Parsons, 2013). The formulation of (Baby, Gole, & Lavergnat, 1988) is then used to calculate the zenith total delay, which can be converted to the slant total delay (STD) that is sampled by the radar wave: where θ is the incidence angle of the satellite at each pixel (Doin et al., 2009). The first term relates to the hydrostatic delay and is calculated using the specific gas constant for dry air (R d ), gravitational acceleration at ground level (g 0 ) and pressure (P) at the elevation of ground level (z 0 ). The second term is integrated between ground level (z 0 ) and a reference height (z re f) of 15 km. This term relates to the wet delay and is calculated using the partial pressure of water vapour (P w ), absolute temperature (T) and the inverse compressibility factor for wet air (Z w − 1 ). We assume Z w − 1 = 1 i.e. negate compressibility (e.g., Baby et al., 1988;Doin et al., 2009). The remaining parameters k 1 , k 2 ′ and k 3 are atmospheric refractivity constants from (Smith & Weintraub, 1953). The resulting vertical profiles of STD are horizontally interpolated to the resolution of a digital elevation model (90 m SRTM: Farr et al., 2007) using a bilinear interpolation. For each pixel the correct altitude is selected to produce a 2D map of atmospheric delay measured in millimetres, hereby referred to as an atmospheric phase screen (APS). In this study we take the simplest and most easily transferable approach and estimate all components of the atmospheric phase screen using atmospheric reanalysis data only. These are the only datasets to provide all of the required fields (temperature, pressure, specific humidity) in a physically consistent way. Other archived data such as MERIS and MODIS also provide global measurements of atmospheric water vapour, albeit only in cloud free, day-time conditions, and therefore need to be combined in a self-consistent manner with measurements of pressure and temperature from other sources to calculate STD using Eq. (1) (Pinel et al., 2011;Walters et al., 2013). In response to increases in the use of large-scale atmospheric models to correct InSAR data, systematic global validations between the wet delay retrieved from ERA-I and MERIS are being undertaken (Walters et al., 2015). These results will aid in estimating the uncertainties on the wet delay calculated using ERA-I, and will therefore help to identify regions where it may be beneficial to use data from other sources.

Case study volcanoes
In this section, we evaluate the use of ERA-I and NARR for correcting atmospheric errors in interferograms. First we compare the performance of ERA-I and NARR using the InSAR dataset covering Lassen Volcanic Center. To correct each interferogram, we produce a differential APS for the interferometric pair using the method in Section 3.2 and then subtract this from the interferogram (Fig. 2A). Long wavelength orbital errors are accounted for by solving for linear ramps either before correction for the uncorrected interferograms, or after correction for the corrected interferogram (Biggs et al., 2007;Gourmelen, Amelung, & Lanari, 2010). We compute the reduction in interferogram standard deviation after removal of each APS to determine which model is most applicable in this setting (e.g., Jolivet et al., 2014;Li, Ding, et al., 2006).
We use the chosen model, and interferograms from Lassen Volcanic Center and Medicine Lake Volcano, to assess how the effect of the model prediction relates to different types of atmospheric artefacts (stratified or turbulent) after Jolivet et al. (2014). To systematically differentiate between the two types of noise, we plot the pixel delay of each uncorrected interferogram against corresponding elevations from a digital elevation model and calculate the correlation coefficient, R 2 (Section 2.2.1). Larger R 2 values (N0.1) indicate that interferograms are dominated by topographically-correlated atmospheric artefacts, and therefore low levels of atmospheric turbulence. We then compare the R 2 values to the reduction in standard deviation (Jolivet et al., 2014), and describe any artefacts in the interferograms that are not corrected by atmospheric models.

Comparison between ERA-I and NARR
Using the dataset covering Lassen Volcanic Center, we compare interferogram standard deviation (σ) prior to and post correction with NARR and ERA-I and identify cases where: a) σ is reduced by both models; b) σ is increased by both models; and c) where σ is reduced by one model and increased by the other (Fig. 2A). We find that σ is reduced in 79% of interferograms corrected by NARR compared to 58% of interferograms corrected using ERA-I (Fig. 2B). The mean σ reduction for NARR is 22% and for ERA-I is 18%. For cases where σ has been increased (21% of interferograms for NARR and 32% for ERA-I), the increase in σ is 8% for NARR compared to 14% for ERA-I. Overall we therefore find that NARR reduces interferogram σ more than ERA-I in this setting (Fig. 2B), most likely due to the higher temporal resolution of data assimilations (3 h compared to 12 h for ERA-I) and higher spatial resolution of model grid nodes (Fig. 1). We therefore select NARR to use for further analysis. These findings are comparable to those of Jolivet et al. (2014) who, using an interferogram at Kilauea volcano, found NARR reduced the standard deviation by 83% compared to 27% for ERA-I.

Correcting stratified and turbulent atmospheric delays using NARR
At Lassen Volcanic Center, we find a correlation between atmospheric stratification and the reduction in σ (Fig. 3A). This agrees with our observation that stratified atmospheric artefacts dominate interferograms at this volcano. These atmospheric conditions are well replicated by the model (see example in Fig. 3B), as they are stable over wavelengths comparable to the model-node spacing (32 km), and the resulting artefacts closely resemble the topography from the digital elevation model. When we observe an increase in σ, this may be due to NARR mismodelling atmospheric processes, particularly the distribution of atmospheric water vapour (Doin et al., 2009).
At Medicine Lake Volcano we find that, although there is an overall reduction in σ for 60% of interferograms, NARR has a relatively minimal effect upon the data, reducing σ by 8% on average and by a maximum of 22% (Fig. 3A). We do not observe a correlation between the extent of atmospheric stratification and σ reduction. The topography of Medicine Lake Volcano is less significant than that of Lassen Volcanic Center, and therefore the extent of topographically-correlated atmospheric artefacts is limited. This agrees with visual inspection of interferograms, and past InSAR studies, that suggest turbulent atmospheric artefacts dominate in this dataset (Parker et al., 2014). Turbulence in the lower troposphere is not predicted by NARR using Eq. (1) (Doin et al., 2009;Jolivet et al., 2011Jolivet et al., , 2014 and we therefore do not expect atmospheric artefacts in data from Medicine Lake Volcano to be removed using this method. The atmospheric artefacts in interferograms from Medicine Lake Volcano have quasi-systematic characteristics (Fig. 3B), which have been referred to as ripples (Li, Fielding, Cross, & Muller, 2006) or atmospheric rolls (Parker et al., 2014) and have a typical wavelength of 4-12 km. Similar features in interferograms across the Los Angeles basin are attributed to gravity waves (Li, Fielding, et al., 2006). Medicine Lake Volcano is located east of the main Cascades axis, and Mount Shasta (elevation~3.5 km), lies between Medicine Lake Volcano and the Pacific Ocean in the path of the prevailing wind (Fig. 1). When wind passes such a topographic obstacle, a train of downwind leewaves may be created with an average wavelength of 2-40 km (Price et al., 2013). This phenomenon is not predicted at the resolution of NARR (~32 km), but by incorporating other parameters such as wind speed, there may be scope to systematically identify conditions that would introduce such atmospheric artefacts.
These results suggest that uniformly correcting all interferograms using NARR may result in spurious results at some volcanoes, and that the effects will depend upon the prevailing atmospheric conditions. For volcanoes where topographically-correlated atmospheric artefacts dominate (e.g., Lassen Volcanic Center), correcting all interferograms may reduce the period of time required to measure deformation at a given rate, whereas if turbulent atmospheric conditions prevail, (e.g., Medicine Lake Volcano) this would not necessarily be the case.

Arc-wide assessment of atmospheric uncertainties
InSAR studies commonly investigate uncertainties due to atmospheric stratification using the gradient between the elevation and delay of interferogram pixels (e.g., Delacourt et al., 1998;Taylor & Peltzer, 2006), and estimate temporal variations in atmospheric delays using time-series methods (e.g., Ebmeier et al., 2013a). By using APS data rather than interferograms, we show how these approaches can be used a priori to independently estimate the magnitude of atmospheric uncertainties on an arc-wide scale. Implementing this for all volcanoes in the Cascades, we assess the influence of topographic and geographic variables upon atmospheric uncertainties by comparing these values to the elevation of each edifice, elevation of surrounding where dates are in yymmdd-yymmdd format. Corrected interferograms are labelled with the standard deviation (σ) reduction and an arrow indicating a decrease or increase. Example a) Interferogram σ is reduced by both models. Example b) Interferogram σ is increased by both models. Example c) Interferogram σ is decreased by NARR only. Black lines show the outline of Lassen Volcanic National Park. B) Comparison between σ reduction for the full set of interferograms that has been corrected using ERA-I and NARR. Dashed line is a 1:1 line -above the line NARR correction is better, below the line ERA-I correction is better. Circles represent interferograms, with red circles indicating the examples used in A. Shaded regions are labelled to highlight that more interferograms have an σ increase after correction with ERA-I than NARR. region, relief of each edifice, surface area of each edifice above a threshold elevation, and latitude.

Method
Elevation-delay gradients and time-series of atmospheric delay are calculated by using NARR to produce 100 km × 100 km APSs centred on each edifice, simulating atmospheric conditions for an ascending SAR acquisition on the first day of each month throughout 2009 and 2010 (Section 3.2).
Elevation-delay gradients are calculated by plotting the simulated atmospheric delay from the APS against elevation from corresponding pixels of a digital elevation model (Section 2.2.1), and using linear regression to solve for the best-fitting straight line. Although the relationship between altitude and atmospheric water vapour is approximately exponential (Foster & Bevis, 2003), linear fits have been shown to provide a good approximation over small scales (e.g., Cavalié et al., 2007;Elliott et al., 2008;Wicks et al., 2002). The elevation-delay gradient at each volcano is the mean gradient calculated over the two-year test period (Fig. 4).
The temporal variability of atmospheric delay at each volcano is calculated using time-series of the delay within 3 km of the volcano summit relative to a reference annulus at 15-20 km distance (Fig. 5A). The RMS of the time-series is used as a measure of the temporal variability of atmospheric delays, and we compare the values from the Cascades to those of Ebmeier et al. (2013a), who applied a similar method to interferograms from the Central American Volcanic Arc (Fig. 5B).
Finally, we use the elevation-phase gradients and time-series RMS to investigate the influence of geographic and topographic variables on atmospheric uncertainties. Topographic variables are calculated using a digital elevation model: the absolute elevation of the edifice is found by averaging the elevation of pixels within 3 km of the volcano summit; the elevation of the surrounding region is calculated by averaging the elevation of pixels within an annulus of 15-20 km from the summit; the relief of the volcano is calculated by differencing the summit and reference elevations; and the volcano "mass" (surface area of the edifice Fig. 3. A) Assessment of the correlation between topographically-correlated atmosphere errors in uncorrected interferograms and reduction in interferogram standard deviation (σ) after correction. NARR significantly improves σ for Lassen Volcanic Center but has little impact upon data from Medicine Lake Volcano where atmospheric turbulence dominates. B) Digital elevation models and example interferograms for each volcano. Example for Lassen Volcanic Center (070712-071025) demonstrates topographically-correlated delays, with Lassen Volcanic National Park shown by the black line and Lassen Peak marked by a triangle. Example for Medicine Lake Volcano (080515-080619) demonstrates both topographically-correlated and turbulent delays, with the extent of lava flows and caldera shown by solid and dashed black lines respectively. The direction of the prevailing wind is shown by the blue arrow as in Fig. 1 . above a threshold elevation) is found using the percentage of pixels located within 15 km of the volcano that have an elevation N80% of the summit elevation (Fig. 6). This threshold elevation is selected to ensure no pixels are within the reference annulus at 15-20 km distance. We also test the use of different radii to calculate the absolute and reference elevations to ensure that the values chosen do not significantly affect the analysis. We plot each topographic and geographic factor against the elevation-delay gradients and temporal variabilities calculated for each Cascade volcano and calculate the correlation coefficient, R 2 , from a linear regression (Fig. 6). A large R 2 value indicates that the variable has a significant influence on atmospheric uncertainties, whereas a small R 2 value indicates that the variable has minimal influence. We discuss any correlations observed, identifying which characteristics are associated with large atmospheric artefacts in the Cascades.

Elevation-delay gradients
The majority of atmospheric water vapour lies below elevations of 2 km (Price et al., 2013), and we observe a significant correlation (R 2 = 0.78) between the elevation-delay gradient and the elevation of the reference annuli (Fig. 6). The correlation with the absolute and relative summit elevations is weaker (R 2 = 0.12 and 0.03 respectively), suggesting that the gradient of atmospheric stratification is greater at volcanoes at lower elevations, with little dependence upon the absolute height of the edifice. The correlation between volcano mass and elevation-delay gradient (R 2 = 0.24) suggests that more massive volcanoes may have smaller gradients of atmospheric stratification, although a larger sample size is required to thoroughly investigate the significance of this correlation. This agrees with the observation that volcanoes with greater landmass at a given elevation will elevate the average temperature (Price et al., 2013), reducing the difference between atmospheric conditions at the base and summit of the edifice.
Elevation-delay gradients vary significantly between the windward and leeward sides of each edifice. At 85% of volcanoes we observe elevation-delay gradients that are up to 15% larger in the windward half of the APSs, where levels of water vapour and precipitation are greater (Fig. 4). The correlation between the elevation-delay gradient and latitude (R 2 = 0.27) shows that volcanoes at higher latitudes (lower temperatures) tend to exhibit steeper gradients of atmospheric stratification (Fig. 6). Thus the effects discussed in Section 2.1 not only affect the coherence of interferograms but also the magnitude of the atmospheric artefacts.

Temporal atmospheric variability
Temporal atmospheric variability is most significantly correlated with the relative elevation (relief) of each edifice (R 2 = 0.71), with weaker correlations for the absolute elevation of the summit (R 2 = 0.64), elevation of the reference annulus (R 2 = 0.10), and volcano mass (R 2 = 0.35) (Figs. 5, 6). Linear regression shows that the gradient of temporal atmospheric variability and volcano relief is 0.7 cm/km. This is lower than a gradient of 1.86 cm/km obtained for volcanoes in Central America (Fig. 5), as equatorial, tropical regions have higher and more variable amounts of atmospheric water vapour (Ebmeier et al., 2013a).  We find that the average temporal variability of atmospheric delays is not linearly related to volcano latitude, which suggests that, whilst volcanoes in the northern Cascades may have larger gradients of stratification, the temporal variability of temperature and water vapour is not dependent upon latitude.

Discussion
Large-scale atmospheric models provide estimates of stratified atmospheric delays in the Cascades, which can be used to reduce topographically-correlated atmospheric artefacts in InSAR data. This is in agreement with studies at the Kunlun Fault, Tibet (Jolivet et al., 2011), Makran, Pakistan and Parkfield, California (Jolivet et al., 2014). Using APS data for the whole Cascades volcanic arc, we have demonstrated that large-scale atmospheric models can also be used to investigate atmospheric uncertainties on regional scales, and we estimate that the magnitude of atmospheric uncertainties in the Cascades ranges between 2.1 cm at Mount Shasta, and 0.16 cm at Crater Lake (Fig. 5, Table 2). Here we show how these independent, a priori estimates of atmospheric uncertainties can be used to define detection thresholds and errors for long-and short-term monitoring strategies on an arcwide basis.

Detection thresholds for InSAR studies at the Cascade volcanoes
InSAR studies commonly use techniques that combine information from many interferograms to reduce atmospheric noise and identify small magnitude deformation signals (see references in Section 2.2).
Using the magnitude of atmospheric uncertainties predicted using NARR, it is possible to estimate the number of interferograms required to detect long-term, small magnitude deformation signals, or the minimum magnitude of deformation detectable during periods of volcanic unrest. Long-term deformation, such as subsidence at Medicine Lake Volcano , may be linear with time, but intereruptive and pre-eruptive displacements are often non-linear (e.g., Three Sisters: Dzurisin et al., 2009) and may be aliased between SAR acquisitions. However, for simplicity here we assume linear rates of deformation to estimate detection thresholds.
Past studies in the Cascades found that poor coherence due to temporal decorrelation limited the use of InSAR (Parker et al., 2014;Poland & Lu, 2008;Poland et al., 2004Poland et al., , 2006Wicks et al., 2002). We therefore focus on multi-temporal methods that use the shortest duration interferograms to minimise the loss of coherence over time. For acquisitions at epochs 1 to N, we would expect to produce N-1 short duration interferograms, where the temporal separation between each epoch (and timespan of each interferogram) is equal to the repeat time of the satellite, t r . Past satellites had repeat intervals of 35-46 days e.g. ERS-1/2, ENVISAT and ALOS, but here we look forward to Sentinel 1A, which will have a repeat interval of 12 days for each satellite, therefore reducing temporal decorrelation. The resulting interferograms may be chain-stacked such that the slave image of the first interferogram is the master of the second (Biggs et al., 2007;Walters et al., 2013). Assuming that the deformation measurements are linearly proportional to the time span of the interferogram, we can then use a weighted least squares problem to solve for a linear deformation rate, dϕ dt : where T = [t r , 2t r , …Nt r ] T are the total duration of the observations between epochs 1 to N, and P =[d 1,2 , d 1,3 , …d 1,N ] T are the InSAR measurements of cumulative displacement acquired over these times, with the subscripts indicating the epochs. The inversion is weighted by Σ p , the variance covariance matrix for the InSAR displacement observations P. The variance of each interferogram, σ i 2 , appears on the diagonal of Σ p . Rather than calculating σ i 2 from a set of interferograms (e.g., Hanssen, 2001), we use the a priori estimates of atmospheric uncertainties calculated using the RMS of time-series of atmospheric delay described in Section 5.3 (Table 2). We assume that each epoch can be treated independently therefore the variance on each interferogram, σ i 2 , is the sum of the variances on each epoch, σ e 2 , such that σ i 2 = 2σ e 2 . As all interferograms share a common master image (the  first epoch), every observation in P has a covariance of σ e 2 with every other observation, and so the off-diagonal terms of Σ p are equal to σ e 2 .
The standard error associated with the linear deformation rate is then, which we can use to estimate operational parameters for long-and short-term volcano monitoring.

Long-term volcano deformation
Long-term, linear deformation rates in the Cascades are measured to be~1 cm/yr. at Medicine Lake Volcano  and Lassen Volcanic Center (Poland et al., 2004). Using Eq. (3), we estimate how many consecutive interferograms would be required to measure deformation at this rate at each Cascade volcano such that σ r b deformation rate. At the repeat interval of Sentinel 1A, between 4 and 11 interferograms are required to measure this small magnitude deformation signal in the Cascades (Table 2) equating to monitoring periods of between 60 days (Crater Lake) and N130 days (Mount Shasta). This is the minimum period of monitoring required to measure such deformation, and a longer observation period would be necessary if interferograms in the chain were not useable due to incoherence. Past studies in the Cascades have found that in some cases only interferograms formed during summer months (July-October) are coherent (e.g., Three Sisters: Dzurisin et al., 2009), which would increase the required monitoring period by up to 1 year (Mount Shasta). The extent to which incoherence compromises the estimates of minimum monitoring periods will vary between geographical regions.

Short-term pre-eruptive unrest
Using InSAR to monitor ground deformation during periods of volcanic unrest avoids exposure of ground-based personnel to the dangers posed by an active volcano (e.g., Dzurisin, 2007). As such, satellite measurements have been used to monitor pre-eruptive behaviour (e.g., Piton de la Fournaise: Peltier et al., 2010) and make decisions during periods of unrest (e.g., Merapi: Pallister et al., 2013).
Time-scales of unrest preceding volcanic eruptions span several orders of magnitude (Passarelli & Brodsky, 2012), and estimating the duration of pre-eruptive unrest at volcanoes requires numerous judgments and assumptions about what characterises unrest. However, studies have identified statistical links between run-up time and composition (Passarelli & Brodsky, 2012) or volcano type (Phillipson, Sobradelo, & Gottsmann, 2013). To define a duration of pre-eruptive unrest for each Cascade volcano we use the volcano type as defined by the Smithsonian Global Volcanism Programme, and a global synthesis of volcanic activity, which shows that unrest preceding eruptions spans a median average of 2 days for complex volcanoes, 1 month for stratovolcanoes, 2 months for calderas and 4 months for shield volcanoes (Phillipson et al., 2013). We compare these values to the repeat time of Sentinel (12 days) to calculate the minimum number of interferograms spanning typical pre-eruptive unrest periods at each volcano type. Using Eq. (3), we then estimate the standard error for chain stacks of these InSAR measurements, and therefore the minimum observable deformation preceding eruption (Table 3). (Note that the unrest periods from Phillipson et al. (2013) are associated with large uncertainties and we include a table in the Supplementary information that reflects the full range of durations.) At stratovolcanoes (which includes 9 of the 13 Cascade volcanoes) we expect per-eruptive unrest to last on the order of 1 month, which would be covered by up to 5 interferograms. Of the Cascade stratovolcanoes, the lowest detection threshold is at Lassen Volcanic Center (b4 cm) and the highest is at Mount Shasta (~7 cm). For calderas, the expected period of unrest is on the order of 2 months, equivalent tõ 5 consecutive interferograms, and equating to a detection threshold of 1-2 cm for calderas in the Cascades. The longest pre-eruptive unrest periods are associated with shield volcanoes (~4 months), which we expect to be covered by~10 Sentinel interferograms. Shield volcanoes therefore yield the smallest detection thresholds of b1 cm. Our analysis suggests that no pre-eruptive unrest would be detectable at complex volcanoes such as Three Sisters Volcano, where unrest periods are on the order of days and therefore less than the repeat interval of any SAR satellites (Phillipson et al., 2013).
The only Cascade volcano to have undergone unrest and eruption in the last century is Mount St. Helens. In the 3-4 weeks prior to the catastrophic May 18th 1980 eruption, ground-based measurements revealed that the "bulge" forming on the northern flank of the volcano deformed at up to 2.5 m/day (Lipman & Mullineaux, 1981). Deformation of this magnitude would result in incoherence, but signal on the periphery of the bulge would likely be measurable above atmospheric noise in a single interferogram. However, when Mount St. Helens reawakened in 2004, deformation rates equated to b1 cm of ground deformation measurable in a single interferogram, and hence InSAR measurements were dominated by atmospheric artefacts (Poland & Lu, 2008).

Applicability to regional volcano InSAR studies
As the temporal coverage of global SAR data increases, large-scale atmospheric models are an accessible method of quantifying and correcting atmospheric errors in InSAR datasets. Almost all global regions are covered by models such as ERA-I, and in many regions higher resolution models (e.g. NARR) are also available. Not only are these model datasets available online, but documented open source software such as PyAPS (Python-based Atmospheric Phase Screen) implements the methods used in this study to calculate APSs (Agram et al., 2013).
Unsurprisingly, our tests with data from the Cascade volcanoes show that local, turbulent atmospheric artefacts are not well replicated by atmospheric models as localised, km-scale features such as leewaves are not reproduced by the coarse model-node spacing. Instead, turbulent phenomena may be simulated using non-hydrostatic models such as NH3D (Miranda & James, 1992), which has been used to compute the forcing of air over and around the orography of Mount Etna (Wadge et al., 2002;Webley et al., 2004). However, this approach not been widely used to correct InSAR data and is best suited to specific targets where atmospheric conditions are well constrained (Wadge et al., 2010) rather than large-scale surveys. Although there remains a need for widely applicable methods of modelling atmospheric turbulence, spatial and temporal averaging utilised by multi-temporal analysis goes some way towards mitigating these phase delays (e.g., Cavalié et al., 2007;Doin et al., 2009;Ferretti, Prati, & Rocca, 2001;Zebker et al., 1997), and it is possible to identify and avoid acquisitions associated with tropospheric turbulence using pair-wise logic (e.g., Massonnet & Feigl, 1995). We suggest that there is also scope to use parameters from large-scale atmospheric models, or elevationdelay correlations as in Section 4.2, to systematically identify cases where turbulent conditions dominate and therefore the corrections used in this study may not be effective. Identifying data in this way would be advantageous as focus shifts towards rapid automated processing and correction of interferograms. At volcanoes across the globe, whether covered by many or few SAR acquisitions, steep topography, extensive vegetation, snow cover, unstable deposits, and the surrounding ocean all contribute to incoherence (e.g., Lu & Freymueller, 1998;Massonnet, Feigl, Vadon, & Rossi, 1996;Rosen, Hensley, Zebker, & Webb, 1996). Where InSAR data are coherent, coverage is often restricted to near field deposits, which may be undergoing deformation associated with deposition (Ebmeier, Biggs, Muller, & Avard, 2014). Poor coherence inhibits the use of empirical estimations of the correlation between elevation and delay, even with the use of advanced spatial filtering techniques. Corrections using external datasets are also limited, as few volcanoes are covered by GPS networks or other local weather data. In these cases, atmospheric models provide a globally applicable alternative to correcting and estimating the magnitude of topographically-correlated atmospheric uncertainties on the scale of volcanic arcs.
InSAR data is being used increasingly in timely evaluations of volcanic hazards (e.g., Pallister et al., 2013), and future developments will involve expanding the use of InSAR from a research tool used after the fact, to an active monitoring tool that can inform hazard management during unfolding volcanic crises (Pinel et al., 2014). Data from Sentinel 1A is expected to be available within minutes of acquisition, but at present, atmospheric reanalysis such as ERA-I and NARR are available only after a 3-6 month delay. Therefore, although large-scale atmospheric models are shown to remove atmospheric artefacts in InSAR data from a variety of settings (Doin et al., 2009;Jolivet et al., 2011Jolivet et al., , 2014Pinel et al., 2011;Walters et al., 2013), using these models to calculate and remove APSs from interferograms in near-real-time during volcanic unrest may be an unrealistic goal. In this study, we present an alternative use of large-scale atmospheric models, and demonstrate how they can be used independently from InSAR data to estimate atmospheric uncertainties and calculate detection thresholds a priori. These estimates of uncertainties can be used to establish errors for time-series analysis, and provide decision-makers with an appropriate representation of the accuracy of InSAR data at a given volcano. Using this approach on an arc-wide scale is also informative for the placement of other geodetic equipment, such as continuous GPS stations, as volcanoes that have higher detection thresholds may benefit more from additional ground-based geodetic equipment than those where uncertainties are predicted to be much lower.

Conclusions
Atmospheric artefacts continue to be the most dominant noise source in InSAR datasets. Here we use large-scale atmospheric models to correct topographically-correlated atmospheric artefacts in interferograms from Lassen Volcanic Center. We find that ERA-I and NARR reduce interferogram σ in 58% and 79% of cases respectively, and that NARR is more suited to reducing stratified atmospheric artefacts due to higher temporal and spatial resolution of the model. Using NARR to assess topographically-correlated atmospheric artefacts throughout the Cascades Volcanic Arc we find that: • Elevation-delay gradients range between 0.09-0.19 cm/km and are largest for volcanoes located at lower elevations (e.g., Mount Baker). At 85% of Cascade volcanoes elevation-delay gradients are also up to 15% larger on the windward side of the edifice. • The temporal variability of stratified atmospheric delays increases by 0.7 cm per kilometre of edifice relief, which is much less than estimates for volcanoes in the Central American Volcanic Arc.
Following this analysis we develop a strategy using large-scale atmospheric models to produce a priori estimates of atmospheric uncertainties on the scale of volcanic arcs. For the application of data from Sentinel 1A to the Cascade volcanoes we suggest that: • A minimum of between 60 and 130 days are required to detect longterm linear deformation at a rate of 1 cm/yr. • During periods of pre-eruptive unrest, ground deformation is most likely to be observed at shield volcanoes (detection threshold b 1 cm) and least likely to be observed at complex volcanoes, where pre-eruptive unrest has historically lasted b1 week. Using this approach we can better integrate InSAR data into long-and short-term monitoring efforts by defining the number of images required for multi-temporal analysis, and establishing the sensitivity of data to topographically-correlated atmospheric noise. Estimating atmospheric uncertainties a priori is particularly valuable for near-realtime monitoring of volcanic unrest, when atmospheric reanalysis data is not immediately available and InSAR data are not coherent enough to use empirical methods alone. With the launch of new dedicated SAR satellites and the development of automated processing regimes, quantifying and correcting for atmospheric uncertainties in this way is a useful step towards using InSAR data in near-real-time.