fluctuation and continuous mass loss over the past two decades (Bevis et al., 2019;Chen et al., 2006;Harig & Simons, 2016;Luthcke et al., 2006;Ramillien et al., 2006). However, the GRACE data has limited spatial (∼300 km in Greenland) and temporal (one-month) resolutions. Other methods include using Global Positioning System (GPS) to track ice movement and crustal uplifting (Bevis et al., 2019;Harig & Simons, 2012;Khan et al., 2014), and building input-output flux models through mass budget methods (Mouginot et al., 2019;Velicogna et al., 2020).

Over the past decades, seismologists have used continuous ambient noise recorded by a pair of seismic sensors to retrieve temporal changes in phase delay times, which are then used to measure near surface velocity changes between two locations (Campillo & Paul, 2003;Lecocq et al., 2014). This technique takes advantage of continuous recording from seismic stations, and was initially applied to monitor changes associated with volcanic edifices and fault zones (Brenguier et al., 2008;Mordret et al., 2010;Sens-Schönfelder & Wegler, 2006;Wegler & Sens-Schönfelder, 2007). Recently, it has been used to estimate relative seismic velocity changes (dv/v) in Greenland. For instance, Mordret et al. (2016) study 2-year seismic records from station pairs deployed in southwest/west Greenland, demonstrating strong seasonal fluctuation in dv/v measurements. They propose a poroelastic model to explain these observations, which may result from surface ice mass loading/unloading and induced strains within the crust. Later, Toyokuni et al. (2018) utilize 4.5-year seismic records from station pairs in different parts of Greenland to monitor both long term trends and seasonal fluctuation of dv/v in different regions of Greenland. However, the cross-correlation measurements from station pairs are typically affected by their low singal-to-noise ratio due to energy scattering through long interdistances between stations, which is critical in Greenland due to its sparse station coverage. Thus, in this study, we collect 11 seismographic stations from different regions of Greenland, and retrieve dv/v by using cross-component auto-correlation of seismic ambient noise. We attempt to study both long term trend (including the "2012-2013 transition") and seasonal fluctuation for different parts of Greenland over the past two decades.

We perform cross-component (north-south and west-east components) auto-correlation for each individual station, except for two stations ANGG and ISOG, which are calculated by using cross-correlation due to their short separation (∼50 km). The advantages of using this auto-correlation technique include: (a) we can specify local dv/v at individual station, instead of a wide sampling region between stations. Since the mass change of the GrIS may have strong spatial variations (Bevis et al., 2019;Mouginot et al., 2019;The IMBIE Team, 2020); (b) we can retrieve correlation functions with high signal to noise ratio, because the seismic phases we measure do not travel long distances between stations. Here, we collect broadband (with the sampling rate of 20 Hz) seismic data for nine stations in Greenland, and long period (with the sampling of 1 Hz) data for stations ILULI and DY2G due to their limited bandwidth. We collect data for all stations with availability ranging from 7 to 21 years (more information for seismic stations can be found in Table S1 in Supporting Information S1).

We use the MSNoise package (Lecocq et al., 2014) to perform noise correlation and measure dv/v. All seismic records are demeaned, detrended and filtered from 0.1 to 1 Hz, except applying 0.1-0.5 Hz bandpass filter for stations ILULI and DY2G due to their limited bandwidth. The next step is to compute auto-correlation functions (ACFs) between two horizontal components. We set the analysis duration as 86,400 s (one day) and cut each daily trace into 1,800 s slices (with a 50% overlap) to perform correlation, which gives us all daily ACFs. We set three times root mean squares (RMS) as extreme limits to suppress outliers (e.g., local seismicity) and a spectral whitening for each correlation slice (1,800 s) is used. The final step is to use the moving-window cross spectrum (MWCS) technique to measure dv/v (Clarke et al., 2011), which estimates the relative time delay (dt/t) between current and reference ACFs in the frequency domain. We set the minimum cross coherence value as 0.75, and measure dv/v by using data in the ACFs above this level, so that we can mitigate temporal variability due to inhomogeneous source distributions. If we assume the change of dv/v is spatially homogeneous, then

We can directly retrieve dt/t by using a weighted linear regression (see more details in Text S1 in Supporting Information S1). Figure S1 in Supporting Information S1 shows an example of seismic noise correlation and dv/v measurement for station NRS. To make a balance between mitigating recording gaps and keeping high temporal resolutions, we test different moving-window stacks for ACFs (with 90, 120, and 150 days in Figures S1c-S1e in Supporting Information S1), and finally choose 150 days stacking for the following analysis. The MWCS filter and other control parameters are listed in Table S2 in Supporting Information S1 and shown in Figure S1b in Supporting Information S1.

To test different frequency response for the dv/v measurements, we first perform a short-time Fourier transform to analyze their time-frequency features. We observe that the strongest signal ranges from 0.1 to 0.4 Hz (Figure S2 in Supporting Information S1). Then we measure dv/v using different frequency ranges (Figure S3 in Supporting Information S1). Larger dv/v amplitudes and higher noise levels can be observed in the lower frequency band (0.1-0.5 Hz), and weaker seasonal fluctuation and long-term trend are observed in the higher frequency band (1-2 Hz). Therefore, we choose 0.1-1 Hz for our following discussions. Furthermore, the selection of measurement windows inside the ACFs is quite important for analyzing and interpreting dv/v measurements (Lecocq et al., 2014).

Here, we apply different windows (dt, time delays) from early to late coda arrivals to test the variability of dv/v measurements. We test three different windows (Figure S4a in Supporting Information S1) on four selected stations, the measured dv/v from testing windows present highly consistent seasonal variation (top panels in Figures S4b-S4e in Supporting Information S1) with our current window (20-70 s). This consistency between different windows suggest that our dv/v measurements mainly reflect subsurface velocity changes, instead of changes from ambient noise source distributions. However, we do not exclude that the most direct (from 0 s) or very late coda (>80 s) portions may bring additional biases to our dv/v measurements. Thus, it is safe to select 20-70 s window portion to achieve robust dv/v measurements for the following analysis. In addition, we also compute auto-correlation sensitivity kernels based on assumptions from Pacheco and Snieder (2005) (More details can be found in Text S2 in Supporting Information S1). This auto-correlation kernel can be described as energy diffusion approximation of multiple scattering wave fields, which only depend on time delays (windows) and location in homogeneous media (Pacheco & Snieder, 2005;Richter et al., 2014). In Figure S5 in Supporting Information S1, the results suggest that our selected window portion (20-70 s) is more sensitive to velocity perturbations at depths shallower than 4 km.

The total GrIS mass balance (TotalMB) mainly includes surface ice mass balance (SMB) and ice discharge (D).

The SMB is balanced by the accumulation and runoffs of the surface ice sheet, while D represents ice discharge and calving at the ice-ocean boundaries. Thus, TotalMB = SMB -D (Bevis et al., 2019;Mouginot et al., 2019).

Here, we use a data set that separately reconstructs SMB and D by using regional climate models provided by the Program for Monitoring of the Greenland Ice Sheet, Denmark (Mankoff et al., 2021). This data set provides results for seven drainage basins (south-west: SW; south-east: SE; central-west: CW; central-east: CE; north-west: NW; north-east: NE; north: NO) in Greenland (Mouginot & Rignot, 2019). For both regional and total GrIS, the seasonal fluctuation of TotalMB is dominated by SMB, while D has very weak seasonal fluctuation even though with large portions in some regions (Figure S6 in Supporting Information S1).

In Figure 1, we observe similar seasonal fluctuation of dv/v for most stations that follows the variation of regional SMB. Here, we perform cross-wavelet transform between dv/v for regional SMB and each individual station, so that we can analyze their time-frequency features as well as the lag times between these two time series (Mao et al., 2019;Torrence & Compo, 1998). Moreover, we use coefficient of variation (CV, the ratio of standard deviation to mean value) to quantify the variability of measured lag time. A smaller CV value suggests a lag time with smaller temporal variability. For example, in the bottom panel of Figure 1c, the strongest signal suggests that the dv/v for station SFJD has a clear one-year cycle, and an averaged lag time of 67 days with respect to the SMB in the SW region. More details about wavelet transform analysis can be found in Text S3 in Supporting Information S1.

Five stations (ILULI, SFJD, DY2G, NUUK, and NRS) in the SW/CW Greenland have averaged lag times of 34, 67, 84, -38, and 35 days with respect to the regional SMB, respectively (Figures 1b-1f). In the NE and CE Greenland, both stations DAG and SCO have the same averaged lag time of 64 days with respect to the regional SMB (Figures 1g and 1h). In addition, in the SE Greenland, the dv/v for station pair ANGG-ISOG has an averaged lag time of 55 days with respect to the SMB (Figure 1i). Table 1 summarizes the lag times between dv/v and regional SMB (with CV values) for these stations in Greenland.

The seasonal fluctuation of the GrIS is mainly due to melting water runoffs in the summertime, and higher precipitation and refreezing in the wintertime (Joughin et al., 2008;Mordret et al., 2016;Mouginot et al., 2019;Rignot et al., 2008). However, this process may not work for the central Greenland, which has much lower temperature during the summertime because of its higher altitude (Figure S7 in Supporting Information S1).

Stations SUMG and ICESG are located in the central Greenland, which cannot be included in any one of the seven drainage basins. Therefore, we perform the cross-wavelet transform between dv/v with local snowfall rates for these two stations, since the snowfall probably dominates the surface mass change in the central Greenland ( Toyokuni et al., 2018). We use the snowfall rate data set from ERA5 (https:// www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5), which is the fifth generation global climate and weather data set provided by the European Centre for Medium-Range Weather Forecasts (Hersbach et al., 2018). Moreover, we note that the dv/v for station NUUK does not have a positive correlation with the SMB (Figure 1e). Thus, we also try to correlate dv/v with the local snowfall rate for station NUUK.

In the central Greenland, the dv/v of stations SUMG and ICESG have averaged lag times of -59 and 82 days with respect to their local snowfall rates (Figures 2b and 2d). The dv/v for station SUMG has responses when there are extreme changes in the snowfall rates, such as 2015-2019, even though its poor seasonal fluctuation (Figure 2b). Station NUUK, which shows a negative correlation with the SMB (Figure 1e), has changed to a normal positive lag time of 44 days with respect to the local snowfall rate (Figure 2c). The better correlation with the local snowfall rate for station NUUK possibly comes from its high urban construction and less ice sheet coverage. It is also possible that coastal water table or other oceanic variables may explain the variations of dv/v, since station NUUK is very close to the coastline.

Table 1 also summarizes the lag times of dv/v with respect to the snowfall rates (with CV values) for these three stations in Greenland.

It has been widely noted that the GrIS encounters extreme warm summers in 2010-2012 that have induced severe mass loss, while the mass loss abruptly slows down due to the cold summer of 2013 (Bevis et al., 2019;Harig & Simons, 2016;Khazendar et al., 2019;The IMBIE Team, 2020). In Figure 3a, we show dv/v time series for stations SFJD, NRS and ILULI in the SW Greenland. Their dv/v present indicative turning around early 2013, which possibly correlates with the "2012-2013 warm-cold transition." To better visualize the trends of dv/v, we use the following equation to fit the dv/v time series:

where t denotes the time with unit in day. The two cosine terms are used to fit the seasonal fluctuation, while the final regression term x 5 t + x 6 is used to represent the linear trend of dv/v. Parameters x 1 to x 6 are estimated 2). Note the increase trend for some SW and central stations after 2013. Gray triangles represent the stations with absence of records.

by using the nonlinear least-squares regression. We use F(t) to separately fit dv/v before and after January of 2013 (black curves in Figure 3a). All these three stations involve the evident negative-positive yearly change rates (dv/ dt) before and after 2013 (with 𝐴𝐴 |𝑥𝑥5| > 4 × 10 -3 in Table 2). In correlation with our dv/v, we calculate the annual ice mass change rates (dM/dt, i.e., the first derivative of the absolute mass change time series) for different regions of Greenland (Figure S8 in Supporting Information S1).

Here, we calculate the annual dM/dt by using a 3-year moving window along the time series of the mass change, which is the same procedure as The IMBIE Team (2020)'s analysis. Figure 3a and Figure S8 in Supporting Information S1 show that the dM/dt of the SW Greenland has the strongest change during 2010-2013, with a peak loss rate of -84 Gt/yr in 2011. Furthermore, we present the summertime (average of June, July and August) North Atlantic Oscillation (NAO) index in Figure 3a, which is based on differences between the subtropical sea-level high pressure and the subpolar low pressure (National Oceanic and Atmospheric Administration, 2022). Strong negative NAO phases indicate above-normal temperatures in Greenland, and vice versa. There are six successive negative summer NAO indexes before 2012 and an abrupt change (with a ΔNAO of +2.3) to positive in 2013, which indicates that Greenland has gone through extreme warm summers until 2012, then entered an abnormal cold summer in 2013. Therefore, the consistent turning trends of dv/v, dM/ dt and NAO index shown in Figure 3a may suggest that the decrease-increase trend of dv/v in the SW Greenland is related to the "2012-2013 warm-cold transition."

In addition, we find a similar turning in the dv/v for station SUMG in the central Greenland (Figure 3b). We also see a single increase trend in the dv/v for station ICESG due to its incomplete records before 2012. For another SW station NUUK, the absence of similar long-term trend (Figure S9a in Supporting Information S1) probably comes from its limited data coverage and better correlation with the snowfall rate instead of SMB. From stations in other regions of Greenland, we just observe relatively weak long-term trends and small values of dv/dt (Figure S9d in Supporting Information S1), which suggest their weak long term mass change rates over the study period. To better visualize the spatial change of the dv/v trend, the stations are color-coded by their dv/dt before and after the transition period in Figures 3c and 3d, respectively.

Here, we just fit the dv/v to the end of 2016. The dv/v in the SW Greenland also shows another sharp drop from 2018 to 2019, especially for stations SFJD and NRS (Figure 3a). This observation agrees with previous studies that the GrIS has experienced another extreme mass loss in the summer of 2019 (Velicogna et al., 2020), with ΔNAO of -2.7 from 2018 to 2019 (Figure 3a). Mordret et al. (2016) measured dv/v between station pairs in the SW/W Greenland, and found dv/v has 2-3 months lag with respect to the GrIS mass change estimated from the GRACE measurements. In this study, we observe stable, continuous and positive lag times of dv/v for most stations with respect to regional/local surface mass changes (Table 1), which is in good agreement with observations from Mordret et al. (2016). For station SUMG, we find an abnormal negative lag time (Figure 2b). This poor correlation probably results from complex local surface mass variation that may bias our measurements. Another possible explanation is that station SUMG is located on a very thick icesheet, and dv/v is less sensitive to deeper processes (see different behaviors of sensitivity kernels in shallower ice layers shown in Figure S5 in Supporting Information S1). For other stations, we note that the lag times of dv/v vary from location to location, especially for stations DY2G and ICESG with large lags of 84 and 82 days, respectively (Table 1). Therefore, in this study, we attempt to explore the implications and causes of various dv/v lag times at different regions.

Previous studies have proposed different processes to explain active crustal deformation, such as poroelastic (pore pressure diffusion) and viscoelastic (viscous flow inside the curst-mantle system) processes (Segall, 2010). Mordret et al. (2016) investigated these two end-member models, and concluded that the poroelastic model is probably more appropriate to explain the seasonal variation and lag time of dv/v due to Greenland's regional geological characteristics (fluids and presence of till layers). This argument comes from the modification of Tsai (2011), and includes a nonlinear relation between seismic wave speeds changes and pore pressure variations due to surface mass loading/unloading. Based on this model, the lag time of dv/v (Δt) against the surface pressure can be attributed to the hydraulic properties of incompetent layer and pore pressure diffusion within the bedrock, which can be expressed as (Mordret et al., 2016):

where z t and K t are the thickness and hydraulic diffusivity of the incompetent layer, which can be considered as deformable subglacial till layer that exists between iceberg and bedrock in glacial areas (Iverson et al., 1997;Truffer et al., 2000). K c is the hydraulic diffusivity of the crust, ω and k are the angular frequency and wavenumber of the surface pressure field. The subglacial till layer is the key to understand ice flow and meltwater inputs (Harper et al., 2017), however, direct observations under ice sheets are always challenging. Following Equation 3, we attempt to estimate z t and K t by using measured Δt through a grid search (more details about the method and parameters can be found in Text S4 in Supporting Information S1). In Figure 4b, we compare the misfit curves for all stations with different Δt. It is notable that stations DY2G and ICESG, which have larger Δt, are separated from other stations. Coincidentally, these two stations are located in the central region and far away from the ice-ocean boundaries (Figure 4a). Tsai (2011) has conducted a synthetic test and suggested that the thickness of the incompetent layer (z t ) tends to predominantly control the lag time of dv/v (Δt). Moreover, there are 32 borehole measurements in the Kangerlussuaq sector (with a maximum 30 km separation and <10 km to station SFJD in Figure 4a) in the SW Greenland, which demonstrate that their sediments over bed is <1 m thick (Figure 4b), even absence for some locations (Harper et al., 2017). While, another field investigation has been taken across the central portion of the NE Greenland (NEGIS in Figure 4a), by using radio-echo sounding, GPS and active-source seismic techniques (Christianson et al., 2014). In contrast, their results demonstrate that a dilatant subglacial till layer (at least 3 m thick, see Figure 4b) spreads over a large portion of the central NE Greenland (Christianson et al., 2014). These two field measurements provide us additional references for our speculation that the larger Δt for stations DY2G and ICESG may come from thicker subglacial till layers in the central Greenland. While, other stations near the ice-ocean boundaries with smaller Δt may indicate weak underlying till due to recent deglaciation and/or variable rates of erosion (Harper et al., 2017). We also test Δt of measured dv/v by using different windows (Figure S4 in Supporting Information S1). For each station, the measured Δt from different windows are basic consistent with each other (Figure S4f in Supporting Information S1). We notice some single early or late windows can lead to biased measurements with large CV values, which is possibly due to different phase types shows two field measurements for regional till layer thickness. The Kan. sector represents 32 boreholes near the Kangerlussuaq sector in the SW Greenland (Harper et al., 2017). The NEGIS represents a combined measurement across the central portion of the NE Greenland (Christianson et al., 2014). In Panel (b), each curve represents grid search misfit values less than 0.5 with the combinations of optimal z t and K t for each station. Station codes, regions and observed dv/v lag times (in Table 1) are labeled. Two black dashed lines denote two field measurements of till layers in Panel (a). More details about the parameters can be found in Text S4 in Supporting Information S1.

or higher noise level. Thus, our current window can largely mitigate single window discrepancies and provide reliable Δt measurements. Here, our interpretation of varied Δt is mainly based on the proposed poroelastic model (Mordret et al., 2016), which is just one possible explanation for our observations at the current stage due to our limited knowledge on local glaciological configuration. We certainly cannot exclude other physical processes that may lead to the seasonal variation of dv/v, such as the changes of subglacial channels or cavities and noise source distribution changes.

Another goal of this study is to explore long-term trends of the dv/v. Mordret et al. (2016) have mentioned that the seismic noise correlation technique may not be sensitive to the long-term trends of ice mass loading. However, our results suggest that this seismic-based technique still likely has potential to monitor the longterm trend of the GrIS, since dv/v presents consistent trend with the first derivative of the ice mass change (dM/dt) during the "2012-2013 warm-cold transition," especially in the SW Greenland. The SE Greenland is another region that is probably affected by this abnormal transition, but there is no corresponding trend in dv/v for station pair AGNN-ISOG. Besides the incomplete data records, another possible reason could be the existence of perennial firn aquifer (PFA), which is a liquid water reservoir that persists throughout the wintertime in the SE Greenland (Forster et al., 2014). If this is the case, the measured mass loss by the regional climate models will be overestimated due to the additional meltwater storage of PFA during warm summers (2010)(2011)(2012). As for two central stations SUMG and ICESG, Toyokuni et al. (2018) performed cross-correlation for this station pair, and obtained a similar increased trend of dv/v from 2011 to 2015. We agree with their analysis that the increase probably comes from the snowstatic pressure accumulation. However, we suggest that the sudden increase trend of dv/v for station SUMG after 2013 is possibly related to the cold transition, which has not been noted by Toyokuni et al. (2018). With lacking of reliable physical mechanisms, we claim that the analysis of dv/v long-term trend is mainly based on our observation of consistency with the mass change rate.

We apply a cross-component auto-correlation technique on continuous seismic records from 11 stations in different regions of Greenland to measure relative seismic velocity variation (dv/v) for each individual station. Our measured dv/v have strong seasonal fluctuation, and have indicative long-term trends in some locations. We investigate the correlation between dv/v with ice mass changes in different regions, and observe that the seasonal fluctuation of dv/v has overall less than 3 months lag with respect to the surface ice mass change or local snowfall rates for most stations. The lag times of dv/v may provide constraints for the thickness of subglacial till layers, and a larger dv/v lag time may indicate a thicker till layer, such as in the central Greenland. In the SW Greenland, the long-term trends of dv/v include an abrupt turning at 2013, which may result from the mass change rate due to the "2012-2013 warm-cold transition," However, we cannot exclude other possibilities that may influence dv/v measurements, such as changes of subglacial channels, cavities and noise source distribution. Our observations demonstrate the potential of using seismic noise auto-correlations to monitor seasonal fluctuation and long-term change of the GrIS.