A century of variation in the dependence of Greenland iceberg calving on ice sheet surface mass balance and regional climate change

Iceberg calving is a major component of the total mass balance of the Greenland ice sheet (GrIS). A century-long record of Greenland icebergs comes from the International Ice Patrol's record of icebergs (I48N) passing latitude 48° N, off Newfoundland. I48N exhibits strong interannual variability, with a significant increase in amplitude over recent decades. In this study, we show, through a combination of nonlinear system identification and coupled ocean–iceberg modelling, that I48N's variability is predominantly caused by fluctuation in GrIS calving discharge rather than open ocean iceberg melting. We also demonstrate that the episodic variation in iceberg discharge is strongly linked to a nonlinear combination of recent changes in the surface mass balance (SMB) of the GrIS and regional atmospheric and oceanic climate variability, on the scale of the previous 1–3 years, with the dominant causal mechanism shifting between glaciological (SMB) and climatic (ocean temperature) over time. We suggest that this is a change in whether glacial run-off or under-ice melting is dominant, respectively. We also suggest that GrIS calving discharge is episodic on at least a regional scale and has recently been increasing significantly, largely as a result of west Greenland sources.

Iceberg calving is a major component of the total mass balance of the Greenland ice sheet (GrIS). A century-long record of Greenland icebergs comes from the International Ice Patrol's record of icebergs (I48N) passing latitude 48 • N, off Newfoundland. I48N exhibits strong interannual variability, with a significant increase in amplitude over recent decades. In this study, we show, through a combination of nonlinear system identification and coupled ocean-iceberg modelling, that I48N's variability is predominantly caused by fluctuation in GrIS calving discharge rather than open ocean iceberg melting. We also demonstrate that the episodic variation in iceberg discharge is strongly linked to a nonlinear combination of recent changes in the surface mass balance (SMB) of the GrIS and regional atmospheric and oceanic climate variability, on the scale of the previous 1-3 years, with the dominant causal mechanism shifting between glaciological (SMB) and climatic (ocean temperature) over time. We suggest that this is a change in whether glacial run-off or under-ice melting is dominant, respectively. We also suggest that GrIS calving discharge is episodic on at least a regional scale and has recently been increasing significantly, largely as a result of west Greenland sources.
2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/3.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
The iceberg calving flux is a major component of the total mass balance of the Greenland ice sheet (GrIS). Past work, using a variety of approaches, has found that one of the components of the GrIS total mass balance, namely the surface mass balance (SMB; figure 1) of the GrIS, has high interannual variability, but with two periods of significant decline over the past century, during 1930-1960 [1] and again in the last decade [1,2]. By contrast, current estimates of one measure of the other key component, the solid ice loss, since the early 1990s, namely the ice flux across the GrIS grounding line, suggest this latter function varies much more smoothly on the ice sheet scale [2] but has increased by approximately 20% since the 1990s. These ice loss estimates exhibit no linear correlation with SMB, except when smoothed and considered over the long term [3]. Three questions are raised by these analyses. Is the iceberg discharge, rather than grounding line ice flux, smoothly varying, given the large background flux, the doubling in estimates of its size in the last decade (cf. [2] with [4][5][6]) and historical [7] and recent palaeoceanographic [8] records of strong variability? How has this discharge varied over longer time scales than can be measured from satellites? What are the major environmental variables driving these changes?
Monthly, counts of iceberg numbers observed south of 48 • N in the western North Atlantic have been compiled by the United States Coast Guard since 1913, absorbing previous reports to the US Hydrographic Service back to 1900. While these numbers have always been acknowledged as estimates [9] and come from a range of platforms evolving over time (ships to aircraft, radar, the use of models and more recently satellites [10]), the series (I48N) is the only long-term estimate of iceberg numbers in the northwest Atlantic available. I48N has a strong, and regular, seasonal cycle (figure 2), with a pronounced peak in numbers in spring to early summer. This is likely to be due to a combination of seasonal peaks in discharge [11], a delay effect from the release of icebergs being trapped in winter sea ice [12] and varying travel paths [13]. However, interannual variability in I48N will be due to a combination of variations in the calving flux, overwhelmingly from Greenland [12], as well as climate fluctuations over the northwest Atlantic. It is known that the former can be very variable from year to year for individual glaciers [8,14,15]. However, the Labrador Sea climate can also change very quickly [16]. Therefore, in this study we explore the long-term variability in iceberg discharge from the GrIS through considering I48N using a combination of ocean-iceberg and nonlinear auto-regressive moving average with exogenous input (NARMAX) system identification modelling both to separate the discharge signal in I48N from the mean ocean-atmospheric interactions following calving and to understand the processes and time scales involved in its variation. The NARMAX approach is described below in the Material and methods, but, briefly, it is a modelling framework which allows the user to construct linear or nonlinear dynamic models between inputs (exogenous variables) and outputs (auto-regressive variables) in the presence of coloured and nonlinear noise.

Material and methods (a) Data
I48N came from the International Ice Patrol's (IIP) iceberg counts over 1900-2008 (http://www. avcen.uscg.gov/pdf/iip/International_Ice_Patrols_Iceberg_Counts_ 1900_to_2011.pdf). They are monthly counts of all separate icebergs larger than a visible length of 5 m observed south of a line extending along 48 • N from the Newfoundland coast to approximately 40 • W (see http://www.ec.gc.ca/glaces-ice/ for eastward extent). By the nature of the observations, they will include icebergs from any origin, but are not able to distinguish the calving origin of individual icebergs. A full discussion of this time series is contained in [9,10,12,17], with discussion of its implications for the Greenland iceberg flux found in [12,18]. It has always been the purpose of the IIP to provide ice hazard information to the major shipping routes into eastern North America [9]. Missing significant icebergs entering these shipping lanes has potentially serious consequences, but since the institution of the IIP there have been no serious incidents rspa.royalsocietypublishing.org Proc. R. Soc.  with icebergs, for vessels heeding IIP warnings [9]. This fact alone speaks to the thoroughness of the surveys leading to the I48N compilation. The Greenland SMB was based on a positive degree day run-off/retention model [1], using a 5 km grid [19]. The monthly North Atlantic Oscillation (NAO; [20]) time series is principally component-based. The Labrador Sea surface temperature (LSST) comes from averaging the Kaplan v2 SST [21], over the Labrador Sea area. The daily forcing of both the ocean-iceberg-seaice model and the SMB model comes from the relevant fields of the twentieth century reanalysis (20CR; [22]). Note, however, that the SMB series is a composite based on the 20CR (1871-1957) and European Centre for Medium-Range Weather Forecasts (ECMWF; ERA40 reanalysis 1958-2001 plus ECMWF operational analysis 2002-2010) temperature, and precipitation datasets that were calibrated/validated against in situ data and spliced together [1]. The composite SMB series agrees very well with independent SMB series for the common overlap period of 1958-2010 (see fig. 9 of [1]).

(b) Ocean-iceberg model
The ocean component of the coupled ocean-iceberg-sea-ice model [23] is a curvilinear general circulation ocean model, with its North Pole displaced to central Greenland to increase horizontal resolution in the North Atlantic and Arctic [24], in this case to less than 0.5 • locally around Greenland, and whose sea-ice component is a thermodynamic model with simple advection. To force the ocean-iceberg model on a daily time step over the period covering the same time span as that of the observed iceberg flux record at 48 • N, we use data from the 20CR. This dataset provides rspa.royalsocietypublishing.org Proc. R. Soc.  estimates of global atmospheric variables at a 3 hourly temporal, and 2 • spatial, resolution. It is derived from surface pressure observations, with observed sea surface temperatures and seaice distributions as reanalysis boundary conditions. We calculated daily averages of the variables required for the ocean-iceberg model's heat, fresh water and wind fluxes at each model gridpoint.
A discussion of the ability of the model to reasonably reproduce mean circulation fields is given in [18,23,24]; the strong variability stemming from daily forcing is illustrated in [25]. The coupled iceberg component is a dynamical and thermo-dynamical iceberg trajectory model [26,27], with a range of iceberg sizes released from the main Northern Hemisphere calving sites [23], carrying a weighting scaled by an estimated iceberg flux. This underestimates the total flux from Greenland by somewhat more than half, because it uses only the main 20 release sites and does not include icebergs from the many small tidewater glaciers, although a range of sources from around the entire Greenland coast are used (see fig. 2 of [23] for both the location and relative sizes of discharges). In addition, estimates of the total GrIS calving or grounding line discharge have changed radically in recent years, from around 200 Gt yr −1 [4,28,29] to over 500 Gt yr −1 [2]. Some of this change is the result of improved methods of estimation, and some is the result of likely recent increases in mean flux.
The iceberg model has been extensively tested previously in the Arctic [26], Antarctic [30] and glacial climatic conditions [31][32][33]. Its behaviour due to variation of the various forcings to which icebergs are subject has been fully explored [27]. The iceberg model has shown itself to be robust in capturing the mean trajectories, limits and densities of icebergs.
The first experiment uses the same flux of icebergs each year, but with a regional variation in release proportions through the year. South of 64 • N, 15% of the annual flux is released in October, 75% in January and 10% in April, to match when sea ice breaks up from southern Greenland fjords (from analysis of Landsat imagery), while north of this latitude the releases are split 15% in April, 75% in July and 10% in October to match peak calving times [11,15]. The second experiment used an annual total flux varying from year to year and scaled according to where D g (t) is the annual flux for a given tidewater glacier, t is time in years, I48N m is the mean observed I48N, D gclim is the climatological mean annual discharge of the glacier [26], t = t for glaciers south of 64 • N with t = t + 1 for glaciers north of 64 • N and c is a scaling parameter (270/96) to reflect the underestimate of the total GrIS calving discharge by the model sampling of release sites. There are few well-documented datasets of iceberg trajectories [13,34,35] but the model is consistent with these (e.g. electronic supplementary material, figures S1 and S2).
(c) Nonlinear auto-regressive moving average with exogenous modelling The NARMAX system identification model [36][37][38] uses a forward regression orthogonal leastsquares algorithm to build models term by term from recorded datasets. This is achieved by using the error reduction ratio (ERR), which shows the contribution that each selected model term makes to the variance of the dependent variable (I48N here) expressed as a percentage, taking account of the noise in the data [38]. The NARMAX method searches through an initial library of model terms, which typically includes linear and nonlinear lagged variables, and selects the most significant terms to include in the final model. NARMAX methods have been widely applied to many problems across a wide range of scientific disciplines, such as engineering, modelling space weather, electroencephalography, visual systems in insects and many others (e.g. [39][40][41]; see ch. 14 of [38] for further examples). Annual data of all independent and dependent variables were used to construct the NARMAX models. Three physical variables that will have an impact on melting processes near a glacier's marine termination are the SMB, representing the glacier's surface run-off/accumulation balance, a measure of ocean temperature (LSST), which affects under-ice melting, and an atmospheric circulation measure (NAO). In all cases, we have chosen large-scale quantities as we are considering mean calving behaviour over extensive areas of Greenland, rather than individual fjords. These three variables can be seen as measures of the three glacier retreat mechanisms recently proposed by Straneo et al. [42]. It is worth noting that previous work has shown that the Newfoundland sea-ice extent is strongly correlated with I48N [12]. In a separate study [43], we show that this correlation is not due to an independent process, but is linked to one or more of the three variables we use here for the NARMAX model.
Initial studies showed that auto-regressive terms did not contribute significantly to the models (that is, terms involving I48N itself); such terms were therefore excluded from the library of terms for the model fits. A wide range of potential lag terms extending over several years were also initially tested, but as the significant contributions consistently came from those under 4 years, this period was used to limit the possible number of terms for computational reasons. Similarly, the order of the polynomials principally contributing to the NARMAX models was explored through seeing which terms arose naturally from the NARMAX procedure. The most significant terms, considering their ERR values, were up to third order, so the final NARMAX model only allowed up to third-order polynomials.
In this work, a procedure with a 30 year sliding window, to allow for a temporally evolving solution, has been developed. A previous analysis using this sort of approach in electroencephalography is described in [41]. We did many initial studies on our dataset which are not discussed in this paper to save space (these are currently being prepared for publication). However, mathematical analysis showed that the information in the data is non-stationary and therefore the full dataset cannot be treated as a homogeneous system. We tried different length data windows to optimize the modelling of the time variation within the data and concluded that a 30 year window was the most appropriate. The window is moved through the data, so that each window's model is created from the data of that particular 30 year window. Time variation will therefore be tracked using this approach. This includes effects that are longer than the window length as these longer term changes will moderate the data through this and neighbouring windows in a way that the sliding window approach can track. sensitivity tests involving combinations of variables, and monthly versus annual data, were carried out. Here, we show only the full annual model, but its results are consistent with all the sensitivity tests. Note that a lag here of 0 years between variables does not eliminate the opportunity of cause followed by effect a few months later.

(a) Non-stochastic nature of I48N
It is first shown that I48N cannot be generated by a purely stochastic process, but should be treated as being forced by external inputs. To support this argument, we generated 100 surrogate datasets, denoted by I48N (s) , from the original raw monthly I48N data by using an amplitude adjusted Fourier-transformed surrogates algorithm [44,45]. For each of the surrogate datasets, a nonlinear auto-regressive moving average (NARMA) model of the form was identified, where f (s) (.) are nonlinear functions corresponding to the s-th surrogate, and e(t) is a noise sequence. The maximum nonlinear degree of the model terms was chosen to be 3. In order to show that I48N is not a purely stochastic process, we calculated the normalized mean squared error (NMSE) for each of the 100 surrogate models. Only three NMSE values out of the 100 surrogate tests are lower than that for the original data ( figure 3). This result statistically shows that I48N cannot be a purely stochastic process, even taking into account possible stochastic errors in the variable, and seeking physical causes for its variation is a valid pursuit. For the rest of the paper, we use the annual averaged I48N for clarity.

(b) Ocean modelling
Whether the variation of I48N is dominated by oceanic or glaciological change can be inferred from seeding the coupled ocean-iceberg-sea-ice model with a constant, or unvarying, net annual iceberg flux from major tidewater glaciers and marine ice shelves around the Arctic. The coupled model then predicts well below the observed iceberg number crossing 48 • N in most years (figure 2) and is poorly correlated with I48N (r = 0.11). Variability in the marine environment through which icebergs travel therefore has little discernible impact on modelled I48N and is not the principal factor leading to the large observed variability in figure 2. This means that variations in the upper ocean dynamics caused by interannual variability in wind, and the air-sea heat and freshwater fluxes, and interannual variation in the melting rate of the iceberg, related to ocean temperature and wave strength, are insufficient to explain the highly variable I48N.
If open ocean circulation and temperature change has not caused the dominant variability in iceberg numbers then there must be some fundamental relationship with iceberg discharge from sources whose icebergs reach 48 • N. It is worth noting that this was the implication of a fjordal model study of iceberg sedimentation as well [46], where the glaciological regime was found to be more important than fjordal sea temperature. To confirm this, another experiment was performed, using the same atmospheric forcing, where the annual calving discharge for each Northern Hemisphere glacier was defined by equation (2.1). The resulting model prediction for I48N (variable discharge line in figure 2) shows very high correlation (r = 0.83). From the oceaniceberg modelling, it is found that over 99% of the modelled icebergs crossing 48 • N originated from southern or western Greenland fjords, or other glaciers within northern Baffin Bay ( figure 4). The ocean model results therefore strongly suggest that variation in Greenland glacial discharge should dominate the underlying fluctuation of I48N. As the largest changes in SMB in recent years have occurred in the ablation zone [47], and therefore around calving glaciers, it is likely that the relatively well-known SMB will be physically related to the poorly known discharge [3]. Previously, it has been shown that there is not a simple linear relationship between SMB and I48N Other' is all of eastern Greenland above 65 • N plus sites in Svalbard and the Russian Arctic islands. [1]. A more complex, nonlinear relationship between this index, SMB, and parameterizations of regional climatic and local oceanic change, taken as the NAO and LSST, was therefore sought using NARMAX modelling.

(c) Nonlinear auto-regressive moving average with exogenous modelling
The results of the full NARMAX model of I48N are shown in figure 2 and demonstrate a high level of fit (r = 0.84; see table 1 for characteristic ERR values, but it should be remembered that these are just two realizations of models from the 79 sliding window intervals). Figure 5 shows     dominant. Thus, in the first half of the twentieth century, annual variation in the late spring peak of I48N was controlled by the GrIS SMB of the previous winter, or that two winters before (see figure 6, SMB panel). As suggested by the example shown in table 1 for the sliding window 1900-1929, the dominant relationship at that time tended towards a simple, linear dependence on SMB. However, in the last few decades, the dominant influences have switched to a mix of nonlinear terms, with lags up to 2 years, but with major contributions from 0 and 1 year lags.

Discussion and conclusion
We hypothesize that the reason for the change in NARMAX model behaviour over the twentieth century, particularly the lag changes, is a combination of changing calving origin of those icebergs reaching 48 • N, declining sea ice off Newfoundland [48] and warming sea surface temperatures in the Labrador Sea and Baffin Bay encouraging calving increases [15,49]. We propose that these factors may have led to a switch during the middle years of the century from an earlier, simpler, situation where icebergs tended to originate from southern Greenland (electronic supplementary material, figure S3), and hence reached 48 • N quickly enough for there to appear to be no lag in annual terms, to one where there is a greater mix of icebergs originating from more distant western and northwestern Greenland (from where modelled travel times to 48 • N are typically 1-3 years; see electronic supplementary material, figure S1), plus some still from southern Greenland. We suggest one reason for this change is the observed decline of about a third in maximum sea-ice area in the Labrador Sea [48] from 1920 to 1970, leading to more icebergs escaping from the sea ice in Baffin Bay, rather than a real change in calving origin. However, recent warming of the LSST (figure 1), and so increased sea temperatures near calving fronts [15,49], is likely to have led to a general increase in average calving rates. In addition, strong local correlation of GrIS SMB and NAO immediately up-glacier from the calving zones (figure 7) supports significant interaction between these variables at this local scale, through atmospheric variability, particularly affecting run-off, which has increased significantly in the last two decades [1], and perhaps basal sliding rates near the glacier front. The rise in importance, in models of I48N, of SMB in recent decades, shown in figures 5 and 6, is also consistent with a climatically driven increase in calving, through more run-off favouring faster outlet flow [50][51][52]. Combining these factors, it is likely that there has been a general increase in average calving rates in recent decades, as suggested by the I48N record of figure 2.
Combining the varied model and data sources used in this paper leads to the conclusion that I48N, even at an annual scale, is strongly related to iceberg calving from Greenland, and that the GrIS's SMB and the LSST are the major causal factors linked to this calving. The implication of this finding is that the GrIS calving discharge has increased quickly, but episodically, in recent decades, just as run-off has [53], and the SMB has been declining [1,2,47]. This is consistent with the episodic nature of retreat records from Jakobshavn Isbrae in west Greenland [7], sudden retreats at Kangerdlugssuaq and Helheim in east Greenland in 2004-2005 [8,14,15], and, more generally, the wide range of retreat rates for glaciers around Greenland during 2000-2010 [54]. From our analysis, it is difficult to separate out the variations over time of the regions which have contributed to this episodic flux. I48N cannot differentiate between sources. However, an indication is given by the coupled ocean-iceberg model variation of the origin of icebergs contributing to I48N. This is demonstrated in the electronic supplementary material, figure S2, which clearly shows the switch from southern dominance to a more mixed modelled source distribution from the 1930s onwards. In the last two decades, the model suggests a change from more important northwest Greenland sources to west Greenland sources. A comparison of marine core ice-rafted debris and modelled iceberg origins in the Denmark Strait, off east Greenland, has shown observational support for inferences from modelled iceberg origins [55]. However, it must be remembered that the imposed calving flux varies in the same way for a given year over the whole of Greenland in the model, so what we see in the electronic supplementary material, figure S3, is a convolution of glaciological and oceanic effects varying over time. Nevertheless, the sustained importance of 1 year lags for all variables in the NARMAX model in the last few  Our results suggest that I48N, suitably scaled, is an effective proxy for the GrIS calving flux, or at least for that from the southern and western quadrants of Greenland (figure 2). The series shows short periods of high flux throughout the last century, with particularly high fluxes during more extended periods centred around 1910, 1940, 1958, 1975 and 1985. The most distinctive part of the series, however, is the period since the early 1990s, with more than a decade of generally sustained high flux, but highly variable from year to year. This result can be compared with the distinctly smooth grounding line flux curves produced by van den Broeke et al. [47] and Rignot et al. [2]. These are likely to match calving discharge over longer timescales, but may be problematic in monitoring year-to-year calving variation, which will most likely be driven by local, calving conditions, as suggested by the wide range of retreat rates around Greenland in the last decade [54]. Decadal means of I48N, scaled as in the right-hand side discharge axis of figure 2, are given in table 2. This suggests that, in the mean, despite the large interannual variability, the GrIS in the first half of the twentieth century experienced similar levels of iceberg discharge leading to I48N. There was then a minimum during the 1950s and 1960s, with an increase from the 1970s onwards, peaking in the 1990s at levels similar to, but exceeding those reconstructed by Rignot et al. [2]. Notably, the most recent period (2000-2008) shows a decrease from these levels by over 50%, but still consistent with calving at an historically high level.
The implicit assumption in table 2 is that east Greenland calving behaves in a similar fashion to that from west Greenland, and possible changes in the mean size of icebergs reaching 48 • N reflect changes in tidewater glacier thickness as well as source region. However, even taking into account observational biases over time, this makes the recent past an unprecedented period for high GrIS iceberg fluxes, consistent with rapid changes in Greenland climate occurring around  this time [1], but with significant interannual variability. This will result in a distinctly variable freshwater input to the northwest Atlantic Ocean, and so a more variable interannually changing effect on its circulation.