Skip to main content Accessibility help


  • Access
  • Open access
  • Cited by 2



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Diurnal seismicity cycle linked to subsurface melting on an ice shelf
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Diurnal seismicity cycle linked to subsurface melting on an ice shelf
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Diurnal seismicity cycle linked to subsurface melting on an ice shelf
        Available formats
Export citation


Seismograms acquired on the McMurdo Ice Shelf, Antarctica, during an Austral summer melt season (November 2016–January 2017) reveal a diurnal cycle of seismicity, consisting of hundreds of thousands of small ice quakes limited to a 6–12 hour period during the evening, in an area where there is substantial subsurface melting. This cycle is explained by thermally induced bending and fracture of a frozen surface superimposed on a subsurface slush/water layer that is supported by solar radiation penetration and absorption. A simple, one-dimensional model of heat transfer driven by observed surface air temperature and shortwave absorption reproduces the presence and absence (as daily weather dictated) of the observed diurnal seismicity cycle. Seismic event statistics comparing event occurrence with amplitude suggest that the events are generated in a fractured medium featuring relatively low stresses, as is consistent with a frozen surface superimposed on subsurface slush. Waveforms of the icequakes are consistent with hydroacoustic phases at frequency $ {\bf \gt} \bf 75\,{\bf Hz}$ and flexural-gravity waves at frequency $ \bf {\bf \lt}25\,{\bf Hz}$ . Our results suggest that seismic observation may prove useful in monitoring subsurface melting in a manner that complements other ground-based methods as well as remote sensing.


A leading explanation for sudden ice-shelf disintegration, as exemplified by the 2002 breakup of the Larsen B, involves the effects of atmospheric warming, increased surface melting and the creation of an impermeable ice-shelf surface that retains large volumes of meltwater which eventually drives hydrofracture (Scambos and others, 2000). According to the theory, surface and englacial warming, extended melt seasons and increased down-slope wind conditions gradually turn the thermally protective firn layer into a bare-ice surface that resists meltwater infiltration (e.g. van den Broeke, 2005; Holland and others, 2011; Trusel and others, 2015; Hubbard and others, 2016; Lenaerts and others, 2017). Build-up of surface meltwater storage capacity (e.g. Glasser and Scambos, 2008; Banwell and MacAyeal, 2015; Kingslake and others, 2017) and interruption of flow features, like rivers and streams, that can remove meltwater volume (Banwell, 2017; Bell and others, 2017), may eventually lead to hydrofracture-driven fragmentation of the ice shelf (e.g. Scambos and others, 2000, 2009; Banwell and others, 2013; Macdonald and others, 2018).

Understanding ice-shelf disintegration, and the ability to anticipate it in the future, requires reliable methods to observe surface and subsurface melting. Unlike the ablation zones in Greenland, where steep surface topography drives immediate movement of meltwater as it is created, flat ice-shelf surfaces allow surface and subsurface meltwater to build-up in bog-like areas where the topography is subtle, but heterogeneous (Banwell and others, 2014). Flat surface topography and the absence of firn also promotes subsurface melting driven by solar radiation that penetrates the surface (Leppäranta and others, 2012; Jarvinen and Leppäranta, 2013; Lenaerts and others, 2017). While the collapse of the Larsen B Ice Shelf in 2002 was preceded by a well-developed surface expression of meltwater (Sergienko and MacAyeal, 2005; Glasser and Scambos, 2008), the collapse of the Wilkins Ice Shelf in 2008 was attributed to the presence of water below the surface that enabled hydrofracture (Scambos and others, 2009). The processes allowing development of subsurface water ponding are thus pertinent to ice-shelf stability.

Over a 3-year period from 2015 to 2017, we conducted a ground-based study of surface meltwater features on the ablating portion of the McMurdo Ice Shelf (McMIS), Antarctica (Fig. 1). One objective of the study was to determine the ice-shelf flexure/fracture response associated with surface water ponding and drainage. A second objective of the study, which we report here, was to determine the seismic signals associated with the water ponding and movement. The heterogeneous nature of the surface topography, albedo and debris cover, coupled with meters-thick subsurface meltwater and slush layers (e.g. Paige, 1968; Klokov and Diemand, 1995), made access to the field area and optimal instrument placement virtually impossible (Fig. 2). Because of these difficulties, our seismometer deployment scheme was modest, restricted to two locations and extended for 63 days at the onset of the 2016/17 melt season (Fig. 1).

Fig. 1. Field setting and location of instrumentation. The seismometer deployed at the ‘dry station’ is in a conventional, firn-covered ice-shelf area of ~50 m thickness. The seismometer deployed at the ‘wet station’ is in a region of bare ice and (some) debris cover displaying complex surface meltwater features including surface and subsurface lakes, ephemeral surface ponds with ephemeral ice cover, streams, and ‘relict’ surface lakes that had open water in previous years. Ice thickness at the wet station is ~25 m.

Fig. 2. Photographs of the McMIS surface melt zone taken in December and January. (a) Surface lake near the wet station (28 January 2017). (b) Bog-like conditions near Rift Tip Lake (17 December 2015). (c) Seismometer deployment at the dry station (23 November 2016). (d) Seismometer deployment at the wet station (20 November 2016).

Previous deployment of seismometers on, or adjacent to, floating ice shelves, has been exclusive to cold, firn-covered accumulation zones where surface temperatures are colder than freezing for most, or all, of the year (Podolskiy and Walter, 2016). These previous deployments have focussed primarily on aspects of ice/ocean interaction that give rise to seismic signals, such as flexural wave propagation (Hatherton and Evison, 1962), ocean gravity-wave forcing (Cathles and others, 2009; Bromirski and others, 2010; Bromirski and Stephen, 2012), rift propagation (Bassis and others, 2007; Banwell and others, 2017), tidally modulated stick-slip earthquakes and slip tremor at grounding lines (e.g. Winberry and others, 2013; Lipovsky and Dunham, 2016), vertical stratification of seismic and hydro-acoustic properties (Zhan and others, 2014; Diez and others, 2016) and tidal flexure (e.g. Hulbe and others, 2016). Where these studies have addressed glacial hydrologic processes, the majority have been focussed on subglacial hydrology (e.g. Walter and others, 2008; Winberry and others, 2009). The present study joins a few other recent cryoseismological efforts (e.g. Chaput and others, 2018; MacAyeal, 2018; Podolskiy and Others, 2018) in attempting to characterize surface conditions such as melting, freezing and changing temperature using seismological methods and data.


The McMIS occupies an embayment surrounded by the mountainous Ross Island Archipelago and is neighbored by the dry-valleys of the Royal Society Range to the West, the Ross Ice Shelf to the South and East, and the McMurdo Sound, into which it flows, to the North (e.g. see Fig. 2 of Dunbar and others (2009)). Both the US McMurdo Station and N.Z. Scott Base are located proximal to our field area. Ice flow of the McMIS is primarily from East to West, traversing a strong micro-climate gradient that extends from the firn-covered Ross Ice Shelf, which rarely has a surface temperature above 0°C, to a zone of strong surface ablation made possible by down-slope winds coming over Mount Discovery and Minna Bluff (e.g. Kuipers Munneke and others, 2018). The Western end of the McMIS is covered with wind-blown terrigenous debris (Glasser and others, 2006) and marine sediment derived from basal freezing and anchor ice deposition (Debenham, 1965; Swithinbank, 1970; Kellogg and Kellogg, 1987). The presence of convoluted medial moraines on the Western McMIS, suggests that the McMIS is largely a relict of a once more widespread ice shelf filling McMurdo Sound (Glasser and others, 2006).

The microclimate gradient and the increasing presence of supraglacial debris (a combination of windblown dust and marine material brought up from the ice-shelf bottom) lead to a major transition from cold, firn-covered ice shelf in the East, to bare-ice/debris-covered surfaces in the West. This transition also involves a change in the surface hydrology: first, a surface that is mottled with puddles and slush-filled depressions that are typically covered by thin frozen ice lids (Paige, 1968; Klokov and Diemand, 1995); second, a region of surface ‘bogs’ and lakes within patches of debris; and third, an area of numerous incised sub-parallel surface streams and rivers (Figs 1 and 2).


The two seismometer stations were deployed ~ 20 km apart on either side of the transition between firn-covered and bare/melting ice-shelf surface, and so we refer to them as the ‘dry station’ and the ‘wet station’ (Fig. 1). Both stations were equipped with broadband instruments (Nanometrics T120), and the wet station was additionally equipped with two 3-component short-period geophones (Sercel L28). The instruments operated autonomously from mid-November 2016 to late January 2017 (Table 1). The dry station's single instrument was deployed in a dry snow pit 30 cm below the surface, adjacent to a heavily travelled over-snow road to Pegasus runway. Three instruments at the wet station were deployed in a 1 m deep hole in the bare-ice surface. Due to uncertainty about how the instruments would behave in a wet, ice-shelf ablation zone and concern about widespread flooding, all three instruments of the wet station were deployed within a 5 m2 area.

Table 1. Instrumentation

* Hydroacoustic noise generated by ship operations in McMurdo Sound dominated the data from 15 January 2017 onward.

The Nanometrics T120 seismometers have an instrument response function that is ideal (‘flat’ Bode diagram with amplitude near unity and phase near zero) for periods from 120 s to frequency 100 Hz (see Data from the T120s was dominated by sub 1 Hz sea swell, buoyancy oscillations (ice-shelf bobbing and rocking), impulsive surface gravity waves from nearby calving and iceberg capsize, and teleseisms from distant earthquakes (Lin, 2017). The L28s have an instrument response function that is ideal from ~ 5 Hz to > 500 Hz (again, featuring a ‘flat’ Bode diagram with amplitude near unity and phase near zero, see The L28 instrument response is also ideal for the observation of icequake activity, which is the main emphasis in this study, because it naturally high-pass filters the ground motion to eliminate the sub 1 Hz signal that dominates the record observed by the T120s. A preliminary display of data, high-pass filtered (passing frequencies above 10 Hz), collected by the two T120s operating at each station is provided in Fig. 3 for an 18-day period. The T120 data are high-pass filtered to eliminate the dominant sea-swell signal at frequency ~ 1 Hz and below, and to reveal the diurnal pattern of seismicity (appearing in the top panel of Fig. 3 representing the wet station) to be discussed below. While the data were easily corrected for instrument response, no elements of the analysis presented here involved signals that were outside of the ideal response range of the sensors (‘flat’ Bode diagram with amplitude near unity and phase near zero). Because the T120s were the only common sensors between the two stations, comparisons between the two stations were done with T120 data. When analyzing the diurnal seismicity cycle (see the following sections) at the wet station, data from the L28s was used because their instrument response naturally suppressed sub 1 Hz sea swell.

Fig. 3. Comparison of ground velocity (μm s−1) seismograms (100 sample per second vertical HHZ channel high pass filtered to > 10 Hz; both collected using T120 broadband seismometers, Table 1) for the wet (top panel) and dry stations (bottom panel) over the 18-day period in late November to early December 2016, when differences between the two stations were most apparent (time measure is UTC). The wet station displays a repeating diurnal pulse of seismicity (~6–12 hr in length), and this forms the focus of interest in the present study. The dry station is much noisier due to a predominance of anthropogenic signals from a nearby over-snow road and airfield. The minimum in signal level from 26 to 27 November (blue span lower panel) corresponds to the Thanksgiving holiday (station recess from work activity) for McMurdo Station, when the anthropogenic noise was reduced. The 8-day period over which the standard deviation of the wet and dry ground velocities are compared in Fig. 5 is indicated in red on the upper panel.

Other instruments (see Table 1) involved in our project included two automatic weather stations (AWS), one deployed for a full-year at a site called Artificial Basin (Fig. 1), and another operated for 3 months at Rift Tip Lake during the period of seismometer operation (Fig. 1). The AWS data provide air temperature and net shortwave radiation measurements at 2 m above the surface for comparison with the seismic data. GPS instruments were also deployed in the area, which recorded tidal elevation data for comparison with the seismic data (see below).


Figure 3 displays high-pass filtered (> 10 Hz) vertical ground velocity seismograms for the wet and dry stations over an 18-day subset (22 November--10 December) of the 63-day total observation window, when the differences between the wet and dry ice-shelf seismicity were most apparent. The wet station's ground motion (top panel, Fig. 3) displays a diurnally repeating, ~ 6–12 hr-long, period of increased seismicity that forms the central focus of interest in this study. This period occurred during the local evening. Following 10 December (not shown in Fig. 3), the diurnal seismicity was intermittent and sometimes absent for days. The pattern of diurnal seismicity is not seen in the dry station's seismogram (bottom panel, Fig. 3).

Inspection of the event waveforms in the two stations (Fig. 4) suggests that the events composing the wet station's diurnal seismic signal are generated primarily by natural sources which have small spatial and temporal scales (judging from < 1 s duration and frequency content of the events). In contrast to the short seconds-scale duration of events at the wet station, the majority of signals observed at the dry station are on the order of tens of minutes in duration. This probably reflects anthropogenic sources related to vehicle traffic on the over-snow road located ~ 100 m from the seismometer vault, and possibly also to activity at the McMurdo Station airstrip on the dry part of the McMIS (Williams Field), ~ 4 km distant. The wet station is separated by ~ 20 km from the dry station and associated anthropogenic activity, so the data show no corresponding anthropogenic signals except in late January when hydroacoustic tremor generated in McMurdo Sound by a US Coast Guard icebreaker and the ships it escorted is observed.

Fig. 4. Example waveforms (100 sample per second HHZ channel, high pass filtered > 10 Hz) for the wet (top panel) and dry (bottom panel) stations. Typical signals at the wet station are sub-second in duration and exhibit frequencies well above 10 Hz. In contrast, the dry station exhibits emergent tremor signals that last for tens of minutes with quasi-symmetric amplitude envelopes suggestive of passing vehicles on the nearby road.

To quantitatively investigate the difference between the dry- and wet-station signals, we computed the standard deviation of each station's vertical ground motion evaluated over 60-second time segments and then used a 20-minute running mean to provide a relatively smooth proxy record for event density and seismicity (we shall refer to this proxy record as ‘seismicity’ from this point on). An alternative measure of seismicity, which identifies and counts individual icequakes and other events is developed below in a discussion of how event number varies with amplitude. Both measures of seismicity were found to identify the diurnal cycle, and the standard deviation proxy is preferable for determining correlations with possible causes.

Figure 5 presents each station's seismicity over an 8-day window (25 November–3 December 2016), chosen as it is representative of the differences between the two stations. Prominent in the seismogram of the wet station is the diurnal cycle identified in Fig. 3. Over the 8-day window shown in Fig. 5, the prominent diurnal peaks are visible on 7 days, but one diurnal peak is conspicuously missing on 30 November. We shall use this missing peak as a crucial test in developing a hypothesis for the cause of the diurnal seismicity.

Fig. 5. Standard deviation of ground velocity over 60-second intervals (μm s−1, 100 sample per second HHZ channel high pass filtered to > 10 Hz) for the wet and dry stations over a representative 8-day period at the start of the melt season (time measure is UTC). The data are then smoothed using a 20-minute running mean filter. The wet station's values have been multiplied by 10 to facilitate comparison with the values from the dry station. The wet station has a prominent diurnally repeating pulse of seismicity which appears each day except on 30 November.

Seismicity of the dry station is about an order of magnitude larger in amplitude than that of the wet station, but lacks the prominent diurnally repeating cycle. Although data on local work activity at McMurdo Station and Scott Base were not collected, the peaks and valleys of the dry station's seismicity are consistent with anthropogenic activity on the nearby roads and airfields. In fact, a prominent window of low seismicity is seen from late 25 to early 27 November 2016 in the dry station's record. McMurdo Station celebrated a Thanksgiving holiday work recess during this period. Efforts to find the short-duration, high-frequency events associated with the wet station's diurnal cycle of seismicity in the dry station's record were unsuccessful.


The diurnal pulsation of seismicity observed at the wet station in November and early December, and intermittently from mid-December to mid-January, constitutes the largest difference between the wet and dry station records. The presence of this phenomenon at the wet station, but not at the dry station, is a simple and strong indicator that it is associated with some aspect of ice-shelf melting and surface hydrology.

Diurnal cycles of seismicity on an ice shelf have been observed before, but not in association with ice-shelf melting and surface hydrology. Hulbe and others (2016) observed diurnal and semidiurnal pulsations of seismicity on the Ross Ice Shelf near the grounding line of the Kamb Ice Stream that correlated with falling tide and was caused by tidal flexure induced fracture (see Fig. 7 in Hulbe and others, 2016).

To test whether the diurnal pulsations at the wet station are caused by the tide, we compared the seismicity with the tide observed at a GPS station located 80 m from the seismometer site. The result is shown in Fig. 6. The periods of strong seismicity begin with the rising tide, and seismicity reaches a maximum during the initial period of falling tide. Thus, the mechanism of seismicity generation applicable in the grounding-line situation described by Hulbe and others (2016) could be at play on the McMIS. The most conspicuous refutation of a tidal cause, however, is the fact that a diurnal seismicity pulse is absent on 30 November 2016, despite the ocean tide being approximately the same as for the days when the diurnal seismicity pulse was observed.

Fig. 6. Comparison of ocean tide elevation (m) and seismicity (μm s−1) at the wet station (see Fig. 5). The absence of diurnal seismicity episode on 30 November 2016 is conspicuous given the relatively strong ocean tide at the time. The wet station's values have been multiplied by 30 to facilitate visual comparison with the tidal elevation. We also note the presence of smaller-amplitude seismicity in the record, typically occurring during the last few hours of each day (UTC time). The focus of our study is on the windows of seismicity of large amplitude which occur during mid-day, i.e., between 11:00 and 18:00 UTC.

We conclude, therefore, that a process associated with ocean tide is not the cause of the diurnal seismicity. We qualify this conclusion by noting that there appear to be other, smaller pulsations of seismicity that occur on a daily basis, most notably during the last 2 hours of the day (UTC time), which could be driven by ocean tide (Fig. 5). These smaller pulses occur just before low tide, so they may reflect some aspect of ice-shelf flexure or other movement (Fig. 6).

In searching for another possible cause of the diurnal seismicity, we found studies which reported diurnal seismicity pulsation below Arctic sea ice, near the floating ice cover of Lake Bonney, Antarctica and in debris-covered alpine glaciers of the Himalaya (Milne, 1972; Carmichael and others, 2012; Podolskiy and others, 2018). In all of these studies, the cause of the seismicity was cited to be plate fracture associated with thermal bending moments on the floating ice cover driven by cooling air temperature. The theory of thermal bending moment induced fracture is described by Evans and Untersteiner (1971) and Bažant (1992).

To test this possibility, we compared the seismicity (described above) with the 2-m air temperature measured at an automatic weather station (AWS, data measured every 15 s, and averaged every 15 minutes) located at Rift Tip Lake, ~3 km away from the wet station (Fig. 1). The result is shown in Fig. 7. The comparison reveals that seismicity pulses seem to correspond with the daily cycle of falling air temperature. Most notable in the comparison is the lack of a seismicity pulse on 30 November coinciding with the only day of the 8-day period in which surface temperature did not fall (i.e., the minimum temperature occurs on the day before, and the temperature rises all day to a maximum occurring the following day).

Fig. 7. Comparison of 2-m air temperature (°C, smoothed with a 4-hour running mean filter) with standard deviation of ground motion at the Wet Station (see Fig. 5). The wet station's standard deviation signal has been multiplied by 500 to facilitate comparison with the temperature record. Diurnal seismicity pulses seem to occur at times of falling air temperature except on 30 November when temperature only increased all day.


To investigate further the relationship between diurnal seismicity and air temperature, we created a simple one-dimensional (1-D) (vertical column) numerical model of heat transfer on the McMIS. The purpose of this model is to determine when thermal bending moments on the ice-shelf surface can induce fracture. Key to this model is its ability to simulate subsurface slush/water layers and the development of a centimeters-thick surface ice lid. Thermal bending moments on this lid, and associated fracture are what we propose to be the origin of the diurnal seismicity cycle at the wet station.

The model computes the time-dependent temperature vs. depth profile, T(z,  t), using surface temperature, T s(t), and sea-water temperature, T sw, as boundary conditions at the surface, z = z s and base, z = z b, of the ice shelf. The ice shelf is taken to be 25 m thick, although it varies within the field area from 10 to 30 m (Rack and others, 2013; Campbell and others, 2017). The model features solar radiative heating, S(z,  t), in the near surface of the ice column. The solar heating is taken from shortwave radiation observations at the Artificial Basin AWS, because its time series is longer than at the Rift Tip Lake AWS, and is assumed to decay with depth below the surface with a skin depth of 0.5 m (Leppäranta and others, 2012; Jarvinen and Leppäranta, 2013). When boundary conditions and solar heating allow, the temperature at the surface and in the subsurface can reach the melting point for ice, which we take to be 0°C. The model accounts for latent heat absorption once T = 0°C by changing the water/ice mixing ratio, Q(z,  t), which is 0 for all ice and 1 for all water. Because of the small ice-shelf thickness, heat transfer associated with advection driven by strain rate, surface ablation (e.g. Glasser and others, 2006) and basal platelet ice accumulation (e.g. Smith and others, 2012) are considered second-order effects, and so, for simplicity, only heat conduction associated with shortwave solar absorption is treated here.

A schematic representation of our simple model domain is shown in Fig. 8. In this simple geometry, a frozen lid covers a slushy, water-covered surface when there is sufficient incoming solar radiation to support melting. This concept is based on observations of supraglacial lakes in Dronning Maud Land, Antarctica, in a flat area near a nunatak (Leppäranta and others, 2012; Jarvinen and Leppäranta, 2013).

Fig. 8. Schematic depiction of 1-D model domain used for simple thermal modeling. The temperature- depth profile T(z, t) (left side of figure) represents conditions around ~1 January when solar absorption below the surface generates a subsurface meltwater/slush layer. The schematic depiction of solar heating S(z, t) (right side of figure) displays exponential decay with depth below the surface with a skin depth, d, of 0.5 m, which is taken from observations by Jarvinen and Leppäranta (2013).

The governing equation for the temperature-depth profile is

(1)$$ {\displaystyle{\partial {\cal W}} \over {\partial t}} = k {\displaystyle{\partial^2 T} \over {\partial z^2}} + S(z,t) $$

where t is time, z is the vertical distance (positive upward, 0 at sea level), k = 2.22 W m−1 K−1 is the thermal conductivity of ice, S(z,  t) is the solar heating, ${\cal W}$ is the total heat density (sum of sensible and latent heat) of the ice defined by

(2)$$ {\cal W}(z,t) = \rho_i \ c \ T(z,t) + \rho_i \ L\ Q(z,t) $$

where T(z,  t) is temperature, 0 ≤ Q(z,  t) ≤ 1 is water mixing ratio, and ρi = 917 kg m−3, c = 2.05 kJ kg−1 K−1, and L = 334 kJ kg−1 are the density, heat capacity and latent heat of fusion of ice, respectively. Disregarding pressure effects, we take the melting temperature to be constant, 0°C. This then means that

(3)$$W(z,t) = \left\{ {\matrix{ {\rho _i c\,T(z,t)} & {{\rm if}\quad T \lt 0} \cr {\rho _i\,L\,Q(z,t)} & \,\,{{\rm if}\quad T = 0.} \cr } } \right.$$

To compute T and Q, we use an implicit time-stepping scheme (which is unconditionally stable) with a time-step dt = 30 minutes consistent with the sampling rates of the AWS system providing model forcing data. The vertical resolution of the finite-difference grid, dz, is 3 cm.

At the end of each time step, grid points where T > 0 are adjusted to restore T to 0 by increasing Q by the amount needed to convert the unphysical sensible heat (ρi c T where T > 0) to latent heat (ρi L Q with T = 0). Grid points where Q < 0 are adjusted to restore Q to 0 by decreasing T from the melting point (0°C) to the new temperature T = LQ/c. The water fraction Q is not allowed to exceed 1. When this happens in the model, Q is restored to 1 (at the end of an offending time step) and the extra heat that caused Q to exceed 1 is discarded, as our model is intentionally too simple to account for heat transfer in liquid water layers. Strictly speaking, the adjustments to Q needed to keep it within the [0,  1] interval mean that the model does not conserve energy. Adjustments to prevent Q > 1 were relatively rare and probably acted to reduce the net meltwater production, thus meltwater generation estimates made by the model represent a lower bound.

Boundary conditions required by Eqn (1) are the specifications of the basal and surface temperatures. For basal temperature at z = z b, T(z b,  t) is constant at T sw = −1.9°C, the melting temperature in saline ocean water (e.g. Smith and others, 2012). For surface temperature at z = z s, we specify T s = T(z s + 2,  t) using a derived annual record of temperature measured 2 m above the ice-shelf surface by the AWS located ~ 12 km away at the Artificial Basin (Table 1, Fig. 1). When T(z s + 2,  t) was >0°C, we substituted T s = 0.

To compile the derived annual record, we noted that the AWS at Artificial Basin observed T(z s + 2,  t) for 365 days over the period of 24 January 2016 to 22 January 2017 (this was not the case for the AWS at Rift Tip Lake, which was closer to the wet station). The AWS observations covered 364 days of a 365-day year; thus we assigned each day of AWS data to the appropriate day-of-year and interpolated for the one missing day to derive the annual record.

The derived annual record of temperature used to force the model, T s = T(z s + 2,  t), is shown in Fig. 9. To account for the temperature difference between the Artificial Basin and the wet station site, located 12 km away, a warming offset, T off = +1°C was added to the 2-m temperature shown in Fig. 9. This offset temperature is the value (rounded to the nearest degree) obtained by comparing the AWS stations at Artificial Basin and Rift Tip Lake. (The AWS at Rift Tip Lake operated only from mid-November 2016 to mid-January 2017, so it could not be used to develop the annual surface-temperature record used by the model.) Sensitivity of the results to this offset temperature is discussed below.

Fig. 9. Temperature (°C) measured 2 m above the ice-shelf surface at Artificial Basin, approximately 18 km from the wet station from 25 January 2016 to 23 January 2017 . This record has been re-expressed as a function of day-of-year for a perpetually repeating 366-day leap year to develop the annual record used as the surface boundary condition for the thermal model. This annual boundary condition record was repeated during model spin-up to eradicate the effects of an arbitrary initial condition. Also shown are time-lines for AWS and seismometer operation relative to the derived annual record.

Solar heating, S(z,  t), required by Eqn (1) is specified as an exponential function of depth in the ice,

(4)$$ S(z,t) = s_{{\rm enh}} \times s_{{\rm net}}(t) {\displaystyle{1} \over {d}} \exp{\left( \displaystyle{z-z_{\rm s}} \over {d}\right)} $$

where s net(t) is a derived annual record of the net shortwave radiation absorbed by the ice shelf developed from the radiation measurements made by the AWS at Artificial Basin. This derivation was made in the same manner as used to derive the annual record of surface temperature described above. Other parameters appearing in Eqn (4) are d = 0.5 m, the shortwave absorption e-folding decay scale (Leppäranta and others, 2012; Jarvinen and Leppäranta, 2013), and s enh = 1.5, a non-dimensional enhancement factor that accounts for the increase of solar radiation absorbed in the wet, debris-covered area of the ice shelf near the wet station relative to the Artificial Basin. As for the offset temperature T off used in the specification of surface temperature, the solar radiation enhancement factor was based on a comparison of shortwave absorption between the AWS measurements at Artificial Basin and Rift Tip Lake. Sensitivity to this choice of s enh is discussed below. The integral of S(z,  t) over the range z ∈ [z s,  z b] is ~ s net because d ≪ H. This means that the function expressed in Eqn (4) effectively deposits the entire shortwave radiation absorbed into the several meters of the ice column just below the ice-shelf surface. The annual climatology of s net(t) derived from the AWS is shown in Fig. 10.

Fig. 10. Absorbed shortwave radiation (W m−2) measured 2 m above the surface at Artificial Basin using a 4-component radiometer. The longwave radiation was observed but not used in the parameterization of solar heating within the ice. This record has been re-expressed as a function of day-of-year for a perpetually repeating 366-day leap year to develop the annual record used as the solar forcing parameter s net for the thermal model. This annual forcing record was repeated during model spin-up to eradicate the effects of an arbitrary initial condition. Also shown are time-lines for AWS and seismometer operation relative to the derived annual record.

Model results

The thermal model was run for a period of 21 years, with each year's forcing being exactly the same time sequence of temperature and solar radiation as shown in Figs. 9 and 10. This model run was sufficient to stabilize T and Q so that they repeated the exact same time evolution in subsequent years. The initial condition was an arbitrary linear function between − 10°C at the surface and − 1.9°C at the bottom.

Results of the model simulation are shown in Fig. 11. The model produces temperate, water-filled ice immediately below a frozen lid, that is maintained by heat loss to the atmosphere above. Notably, from mid-December through late January, the surface temperature is at 0°C, and the lid thickness is 0. The net-water thickness (ice equivalent), distributed in the ice/water mixture reaches a maximum of 161 cm. Our simple model does not allow this water to move horizontally, for example, to drain away via surface streams or through the ice-shelf base by moulins. The water remains in place and is thus available to supply latent heat to the ice column through the Austral winter period.

Fig. 11. Results of the baseline thermal model, total water in the ice column (top panel), thickness of frozen lid (middle panel), bending moment on frozen lid (bottom panel). Period when no subsurface (or surface) water exists on the ice shelf is indicated with gray shading. On day 61 (1 March 2016, vertical line in top panel) of the year, the rift on the McMIS that had previously terminated at Rift-Tip Lake (Fig. 1) propagated another 3 km. Inspection of the rift extension on approximately day 315 of the year ( ~11 November 2016) prior to the onset of surface or subsurface melting revealed that, over the Austral winter (probably immediately following the rifting event), subsurface water drained through the exposed walls of the rift, to form icicles seen in Fig. 12. The 8-day period shown in Figs 6, 7 and 13 is indicated on the bottom panel (i.e., days 330–337).

Although there are no direct observations to verify the model results given for the presence of subsurface meltwater, there is indirect evidence available from Banwell and others (2017). They observed a rift propagation event occurring on 1 March 2016 (day 61 of the year) by satellite, and later went to inspect the newly extended rift on 10 November 2016 (day 315 of 2016), which was prior to the onset of the melt season following the rifting event (see timings of rift event and inspection indicated in Fig. 11). As shown in Fig. 12, icicles were found protruding from the rift walls just below the surface of the ice shelf on 10 November. Given that the rift walls were made on day 61 of the year, during the part of the year when the surface air temperature was falling, it is likely that these icicles represent residual subsurface water that leaked out of the rift immediately after its formation. The top panel of Fig. 11 shows that the total water content of the ice column was still above ~ 60 cm, and the middle panel of Fig. 11 shows that the subsurface water was covered by ~30 cm of frozen ice. These results are consistent with what is seen in the photograph of the rift walls shown in Fig. 12.

Fig. 12. Icicles protruding from subsurface layer along walls of a recently propagated rift (opening on 1 March 2016, see Banwell and others (2017)). This photograph was taken on 10 November 2016, which was before renewed surface or subsurface melting on the ice shelf, immediately after the ensuing winter. We conclude that the most likely time for the icicles to have developed was immediately after the splitting of the ice shelf to form the rift walls on 1 March 2016. This was at a time when surface temperatures would have prohibited melting, hence the source of the water must have been subsurface residual from the preceding December/January period.

Thermal bending moment on frozen lid

The ‘frozen lid’ produced by the model is taken to be the cold ice layer with zero water content (Q = 0) that overlies a layer of temperate ice that contains water within its inter-grain space (Q > 0). If there is no temperate/slushy ice below the frozen lid, the frozen lid is assumed not to exist. The frozen lid first appears in the melt season around mid November (~ day 320) as a thin, ~ 0.5-m thick, superimposed ice cover that fluctuates in thickness until early February (~ day 40, Fig. 11). After about day 40, the lid monotonically thickens to ~125 cm as the subsurface water slowly freezes with the end of the austral summer melt season. The lid ceases to exist once the subsurface layer of temperate/slushy ice disappears around day 100 (9 April).

To pursue the hypothesis that the seismicity of the wet station is associated with thermal fracture, we compute the thermal bending moment per meter width of the frozen lid, M T (N m)/m, using the expression developed by Baz̆ant (1992; Eqn 5):

(5)$$ M_T(t) = {\displaystyle{E} \over {1 - \nu}} \alpha \Delta T(t) h_{{\rm lid}}(t)^2 I_T $$

where E = 5 GPa is the Young's modulus for ice, ν = 1/3 is the Poisson ratio for ice, α = 50 × 10−6 K−1 is the thermal expansion coefficient, h lid(t) is the frozen lid thickness (m), I T = 1/12 is a nondimensional shape function, and ΔT is the temperature drop between the lower side of the frozen lid, where T = 0, and the ice-shelf surface where T = T(z s,  t) is expressed as a model boundary condition. The formula expressed in Eqn (5) is for a 2-D plate where all bending is confined to 1 dimension. The result of the baseline simulation is shown in the bottom panel of Fig. 11.

The evolution of the bending moment M T(t) over the melt season essentially mimics the frozen lid thickness h lid(t), due to its dependence on the square of the lid thickness (Eqn 5). What is most notable in the bending moment's evolution is the fact that it varies on a diurnal timescale, and can be compared with the seismicity for the period of 25 November--2 December (days 330–337), which is the time when the diurnal pulsation of seismicity contains the conspicuous ‘missing day's cycle’. A close-up of the baseline model results for 25 November--2 December (days 330–337) is shown in Fig. 13. The seismicity is plotted as a function of dM t/dt using the data and model results for this same period in Fig. 14. Notable in the comparison is that seismicity tracks increases of bending moment. This correspondence explains the conspicuous lack of seismicity on 30 November because the bending moment on that day rarely increases (except for a few short windows on the order of a few minutes). We conclude that the bending moment on a frozen lid produced by the thermal model constitutes a plausible explanation for the cause of the seismicity cycle seen in our data.

Fig. 13. Comparison of thermal bending moment on the simulated frozen lid, M T(t), with the standard deviation of seismicity for November 25--2 December 2016. Seismicity tends to correlate with periods of sustained daily increase of thermal bending moment. When the bending moment is constant, falling, such as is conspicuous on 30 November, or undergoing short-time-period changes, there is little-associated seismicity.

Fig. 14. Seismicity as a function of d/dtM T. Seismicity preferentially occurs when d/dtM T.

Reduction of diurnal seismicity during late December

Having used the 8-day window of observations from 25 November through the end of 2 December to develop a hypothesis for the cause of diurnal seismicity, we now examine the entire seismicity record (22 November 2016--14 January 2017) at the wet station to see if the hypothesis successfully predicts the seismicity over a longer time window. (While our seismometer deployment ran through to 20 January 2017, we discard the window from 15 January onward because of hydroacoustic noise in the seismic data associated with icebreaker and ship operations in McMurdo Sound). As shown in Fig. 15, the magnitude of the thermal bending moment on the frozen lid from 11 December 2016 to 14 January 2017 is much less than that before 11 December. This reduction in magnitude is seen to correspond with the elimination of the diurnal pulsations of seismicity, with the exception of a few isolated days. (Notable days in which the diurnal seismicity returned, ~19 December 2016 and 2 January 2017, were during periods when the weather was relatively cold.) This elimination corresponds with the time when the thickness of the frozen lid was generally zero in the model run (Fig. 11). The zero thickness was, in turn, due to the continued warming of the surface as the melt-season progressed. We regard the correspondence between the seismicity and thermal bending moment as further confirmation that periods of increasing bending moment are what cause the diurnal pulsations of seismicity.

Fig. 15. Comparison of thermal bending moment on the simulated frozen lid, M T(t), with the standard deviation of seismicity at the wet station for 22 November 2016--14 January 2017. Diurnal seismicity pulses are absent during the time periods (e.g. starting on 11 December) when the thermal bending moment on the simulated frozen lid became very small or zero. The seismicity is not shown from 15 to 20 January because it was dominated by icebreaker and ship hydroacoustic noise from McMurdo Sound.

Our automated camera system (Fig. 1, Table 1) produced a record of surface conditions at Rift Tip Lake (Fig. 16). Figure 16 shows conditions on 7 December and 13 December, 2016, the 2 days that bracket the period of transition from robust diurnal seismicity to a complete absence of diurnal seismicity. The main difference visible in the two images is the replacement of a frozen surface with ice-free meltwater. This confirms that the frozen lid disappeared at roughly the time predicted by the model, and that this accounts for the fact that diurnal seismicity only happened sporadically after 13 December.

Fig. 16. Automatic camera photographs from 7 (top panel) to 13 (lower panel) December taken at 4:31 local time (UTC + 13). View is of a portion of Rift Tip Lake taken toward the South (Mount Discovery in the background) with an AWS in the middle ground. Conditions shown here are similar to those at the wet station. Prior to 13 December, and most strongly on 7 December, diurnally pulsed seismicity was observed 3 km to the West at the wet station. As the top picture confirms, this was during a period of time when the surface was covered by a frozen lid. After 12 December, the diurnal pulses of seismicity disappeared for days at a time. As the bottom picture shows, this corresponds to when the frozen lid had been melted through.

Dependence of model results on assumptions

We return to the question of how the comparison between diurnal seismicity and model-derived thermal bending moments on a frozen lid depends on the assumed T off = +1°C and s enh = 1.5. In model sensitivity simulations, we evaluated the two cases: (a) T off = +1°C and s enh = 1.0, and (b) T off = 0°C and s enh = 1.5. In case (a), a frozen lid covering subsurface water did not develop until day 333 corresponding to 28 November 2016, so the comparison between bending moment on the lid and seismicity was unsatisfactory. In case (b), without an offset warming, the period of time from 25 November through 3 December 2016, did not have a continuously present frozen lid, and, as with case (a), the comparison between the model-derived thermal bending moment and the observed seismicity was unsatisfactory. We conclude that, while our model and its forcing have been simplified to the utmost, its results are satisfactory for the simple purpose of identifying the cause of the diurnal seismicity.


To characterize the events of the diurnal seismicity cycle in terms of seismic phases is challenging because we have only one seismograph observing them. We examined a large number of individual events selected randomly throughout the window of observation and found that a few general descriptive comments can be made about the waveforms and phases comprising the vast catalog of events observed. Our comments must necessarily be viewed as preliminary until further study can acquire data from more than one seismograph.

In Fig. 17, we show the spectrogram for a 10-second interval surrounding the 20th most prominent event (in terms of vertical ground motion) observed on 8 December 2016 at ~12:34 UTC. At 5 s in the spectrogram's horizontal axis, the 20th event is seen to be a sharply defined vertical swath of energy that is greatest in the 100--250 Hz frequency band, but which extends down to 1 Hz (the lower bound of sensitivity for the L28 geophone sensor, which was recording at 500 sps). This event's energy is consistent with a body wave propagating to the seismometer site from a nearby location.

Fig. 17. Spectrogram of vertical ground motion (log velocity squared, relative color scale) observed with the wet station L28 geophone on 8 December 2016, during a prominent part of the diurnal seismicity cycle. The 20th most prominent amplitude event (event 20, out of > 1.6 × 105 events detected on that day) is shown along with about 25 other small events. Event 20 shows energy distributed across all frequencies, consistent with an undispersed body wave propagating from a source location near the seismometer. The other events all exhibit a low-frequency cutoff that appears to range between ~ 25 and ~ 100 Hz. Below about 30 Hz, there are other signals, detailed in Fig. 18. Horizontal striping is likely an artifact of the MatLab routine (spectrogram()) used to make the spectrogram.

In the 5-second windows before and after the 20th event in Fig. 17, there are numerous smaller events (their ranking among the ~ 150 000 events observed during 8 December 2016 was determined but not of interest). These events show a curious low-frequency cutoff somewhere between 25 and 100 Hz (depending on event), and there are also discrete bands of energy at the low-frequency portion of the spectrogram, below ~30 Hz, which have an apparent high-frequency cutoff. We speculate that the low-frequency cutoff is likely due to the vertical arrangement of ice over water over solid earth rendering the water column a ‘slow’ layer sandwiched between two ‘fast’ layers in terms of body-wave phase velocity (MacAyeal and others, 2015a,b). This arrangement is analogous to the ‘SOFAR’ sound channel in the ocean (Ewing and Worzel, 1948). The events that show the low-frequency cutoff in Fig. 17 are typical of the 10s of thousands of events of moderate to small magnitude, and this suggests that they are propagating as hydroacoustic phases from event sources distant on the order of hundreds of meters from the seismometer. The cutoff frequency, f cut, is related to the vertical thickness of the water layer,

(6)$$ f_{{\rm cut}} = {\displaystyle{c_{\rm s}} \over {d}} $$

where c s ≈ 1500 ms−1 is the speed of sound in seawater near 0°C, and d is the water-column thickness. This relation, and the 25–100 Hz cutoff frequency range observed in the data suggest that d ranges from ~ 15 to ~ 60 m, which is consistent with the bathymetry and ice thickness of the McMurdo Sound and McMIS, respectively.

At low frequency, below 30 Hz, the spectrogram in Fig. 17 displays signals that appear related to the high-frequency events, but which are restricted to frequencies below 30 Hz. A spectrogram of the low-pass filtered vertical ground motion (< 30 Hz) shown in Fig. 18 reveals that the signals tend to be reverse dispersed, with higher frequency components arriving before low-frequency ones. This is also evident in the seismogram of the signals, as displayed in Fig. 19. Phases displaying reverse dispersion are a prominent component of floating ice-sheet seismology, and are observed commonly in lake ice, sea ice and ice shelves (e.g. Ewing and Crary, 1934; Press and Ewing, 1951; Hunkins, 1960; Hatherton and Evison, 1962; Xie and Farmer, 1994; Dudko, 1999, and references therein). A parsimonious explanation of the reverse-dispersed phases seen in the diurnal seismicity cycle is that they are generated when the superimposed ice lid imparts a vertical impulse to the ice-shelf immediately below as it undergoes thermal fracture (e.g. Yang and Yates, 1995).

Fig. 18. Spectrogram of low-pass filtered (30 Hz) vertical ground motion (log velocity squared, relative color scale) observed with the wet station L28 geophone on 8 December 2016, during a prominent part of the diurnal seismicity cycle. Events near 2 and 4 s exhibit reverse dispersion, where high-frequency energy arrives before low-frequency energy. Waveforms of the low-pass filtered vertical ground motion shown here are shown in the seismogram of Fig. 19. Horizontal striping is likely an artifact of the MatLab routine (spectrogram()) used to make the spectrogram.

Fig. 19. Seismogram of the low-pass filtered (30 Hz) vertical ground velocity (μm s−1) at the wet station associated with the spectrogram of Fig. 18. Events near 12:34:01 and 12:34:03 clearly display reverse dispersion.


A notable aspect of the diurnal seismicity cycle is the fact that it involves very large numbers of relatively small events. To provide some sense of what number of events are involved and how their amplitudes are distributed, we attempted to count them. Following Carmichael and others (2012, 2017), we experimented with using a short-term/long-term signal variance detector to determine the onset and amplitude of individual events. Difficulty in the implementation of this detector due to the variability of intra- and inter-event timescales led us to use the MatLab™ function findpeaks() as a means of identifying and measuring the amplitude of single events. After experimentation with the parameters that characterize the event detection algorithm findpeaks(), we selected a minimum peak threshold (1 μm s−1 in the vertical ground motion) and a minimum peak time separation (0.1 s) to exclude multiple peak detections within a single apparent event. This approach thus means that multiple events which occur within a time window of 0.1 s or less will contribute one detection (the first one to occur) to the event statistics. Furthermore, the detection algorithm does not exclude noise features that are unrelated to the diurnal seismicity, and it is difficult to assess the uncertainty in the event statistics that result. Assuming the noise features to be of uniform amplitude during the entire day, the minimum peak threshold (1 μm s−1) in the vertical ground motion chosen for the detection algorithm was higher than the noise features observed during quiescent periods of the diurnal cycle.

Using the 500-Hz vertical ground velocity channel of one of the two short-period instruments deployed at the wet station (the other instrument did not stay level through the observation period), we identified > 2 × 106 events during a 53-day window when the presence or absence of a diurnal cycle could be attributed to natural causes, and not to interference from anthropogenic icebreaker and ship noise. The 500-Hz vertical ground velocity signal was first high-pass filtered with 100 Hz as the pass frequency to emphasize the dominant frequency band of the majority of events determined from Fig. 17. A summary of the event timing and their amplitude is shown in Fig. 20.

Fig. 20. Individual seismic events detected from 100 Hz high-pass filtered vertical ground motion at the wet station during a 53-day period (23 November 2016--14 January 2017) used to compute event amplitude distribution N. Data from 14 to 20 January 2017 were excluded due to the presence of anthropogenic tremor from icebreaker and ship activity in McMurdo Sound. The event detection algorithm was simple, involving peak detection (findpeaks() from MATLAB) with a specified low-amplitude cutoff and time window following a given event detection during which no further events would be detected. The algorithm likely leads to a bias in detecting noise features at low amplitudes that are unrelated to the diurnal seismicity cycle.

Using the event detection algorithm described above, we investigated how the number of events varied with event-amplitude measured at the seismometer site. This approach constitutes a variant to classical frequency-size investigations, which have been carried out routinely since Gutenberg and Richter (1954) noted that, in any population of earthquakes, the logarithm (base 10) of the number N(M) of events with magnitude at least equal to M varies as

(7)$$ \log_{10} N(M) = C - b \ M $$

Because of the mathematical properties of the logarithm function, Eqn (7) also applies with the same b-value to the opposite of its derivative n(M) = −dN/dM, i.e., to the number of events falling into a magnitude bin of constant width centered on a magnitude value M:

(8)$$ \log_{10} n(M) = c - b \ M $$

In practice, the ‘Gutenberg and Richter law’, hereafter ‘GR’, applies remarkably well with b ≈ 1 in Eqns (7) or (8) to most populations of shallow earthquakes. This property was later justified based on scaling laws, for fault zones with a fractal dimension of 2 (Rundle, 1989). A departure from GR, and, in particular, a b-value >1, illustrates a different mechanism of seismogenesis, for example, involving a 3-D medium as opposed to a fault, as found for certain classes of deep earthquakes (Okal and Kirby, 1995), or in the highly brecciated material where a large number of fractures hinders the growth of the source for large seismic events (e.g. McNutt, 2002). For example, a high b-value is generally a proxy for volcanic earthquakes releasing the stress transferred into the volcanic edifice by the ascent of magma (e.g. Yuan and others, 1984).

In the present study, we can only quantify the events through their recorded amplitude A, rather than through a magnitude M, which would have to be obtained by applying a distance correction to log 10A; with only one recording station, we cannot infer epicentral distance and carry out such a correction. Faced with similar circumstances, many scientists have simply fit a GR law to the distribution of recorded amplitudes (using a logarithmic scale)

(9)$$ \log_{10} n(M) = C_1 - b_1 \ \log_{10} A $$

and interpreted their results under the assumption that the observed opposite slope, b 1, is indeed the physical b-value of the original distribution of seismic sources; this was the case, for example, of the first studies, using a single ocean-bottom seismometer, of the micro-seismicity associated with earthquakes along oceanic rifts (Francis and others, 1977). This assumption would obviously be correct if all events were at the same distance from the receiver, and hence their amplitudes identically affected by propagation. We show in Appendix A that it does indeed hold (at least for the smaller earthquakes in the dataset considered) for a number of geometries where epicentral distance varies, but the spatial density of sources is constant.

In this context, our observable ‘amplitude’ A recorded on the short-period seismometer can be regarded as ground velocity, measured in μm/s−1 (L28, Table 1). We proceed to evaluate the population of events as a function of amplitude by binning the data into intervals αi ≤ A ≤ αi+1, with $\alpha _i = \alpha _0 + i \ \Delta \alpha $, and define a cumulative population N as a function of amplitude A as the number of events with an amplitude at least equal to A; note that this definition is independent of the details of the binning procedure. We then use log 10A as a substitute for observed magnitude M and study the distribution of log 10N as a function of M.

Using the entire 53-day observation window to compute N, the result is shown in Fig. 21. The main feature of this distribution is the increasing negative slope of the distribution with increasing log 10(A). While it is tempting to compare the shape of this distribution, with the slope corresponding to minus the ‘b-value’ of the Gutenberg--Richter distribution for earthquakes, this should not be done, because log 10(A), the argument of N, is not corrected for distance, as is the case for M, earthquake magnitude.

Fig. 21. Distribution of N for high-frequency events (>100 Hz) observed at the wet station during the 53-day window (23 November 2016--14 January 2017) when diurnal seismicity was prominent. The observations from 14 to 20 January were excluded due to the presence of ship tremor.

The population N(A) shown in Fig. 21 is characterized by a GR distribution with a slope b 1 ≈ 1.2 for low amplitudes (A < 60 μm/s−1), beyond which an elbow is present and b 1 increases significantly, to values of the order of 3. As discussed in the Appendix, the low-magnitude end of the population implies a similar GR distribution with an intrinsic b 0 = 1.2 for the events recorded. The experiments reported in Appendix A would suggest that the position of the elbow is related to the maximum (true) magnitude m max of the events, to the dominant frequency of the wavetrains (~100 Hz), to the unknown epicentral distance, and to the seemingly unknown anelastic attenuation along the path to the station.

Irrespective of the fractal dimension of the seismogenic material, many factors can affect the b-value of a population of earthquakes. In very general terms, a b-value slightly larger than 1 could be expected in a fractured medium (Mogi, 1962) where faults will have difficulty growing coherently to large dimensions, and for low stresses, as b-values have been found to decrease significantly with increasing ambient stress (Scholz, 1968). Such conditions seem legitimate in the particular environment of the ice shelf. While these preliminary remarks are encouraging, it would be futile to attempt deeper interpretation in the absence of an adequate knowledge of fundamental parameters such as epicentral distances.

A remarkable property of the low-frequency signals described in the previous section (see Figs 18 and 19) is that they obey different population statistics from those shown for the high-frequency signals. (Restricting the event counting routine to low-frequency signals was accomplished by applying the detection algorithm to a low-pass filtered time series.) As shown on Fig. 22, they exhibit a value b 1 ≈ 3 much higher than their high-frequency counterparts, with a much less pronounced, if at all defined, elbow. Their lower dominant frequencies would argue for longer wavelengths, and hence a generally larger source, even though the range of their ground velocities recorded at the receiver is not significantly different from that for the high-frequency events. In this context, we recall that the scaling of strong motion ground velocities (or a fortiori accelerations) with earthquake size in the near field can deviate from simple power laws predicted by source similitude arguments (Anderson and Lei, 1994; Abrahamson and Silva, 1997). In the general framework of a highly fractured source medium, we would certainly expect a decrease of the population of large events, and thus an increase in b-values at the source. As for the weakening or absence of the elbow, it is probably a direct consequence of the higher b-values, as evidenced on Fig. 23: when fewer large events are present the effect of the upper bound m max on the undersampling of large observed amplitudes will be weaker.

Fig. 22. Distribution of N for low-frequency (< 30 Hz), reverse-dispersed events recorded during a 53-day period (23 November 2016--14 January 2017) when diurnal seismicity was prominent. The observations from 14 to 20 January were excluded due to the presence of ship tremor.

Fig. 23. Populations of observed amplitudes interpreted as ‘magnitudes M’ for a dataset of earthquake sources with magnitudes m obeying a GR law with an initial b-value of 0.5 (Left), 1 (Center) and 2 (Right). The earthquakes are distributed with a constant density μ along a linear segment of length L = 10 km, with the receiver located at one end of the segment. The plus symbols are the populations of the individual magnitude bins, and the red dots are the cumulative values. The apparent b-values, b 1, are obtained by linearly regressing the cumulative populations.


One signal revealed by our modest seismometer deployment was not previously observed on an ice shelf and it was only observed in the ‘wet’ area of the McMIS (Fig. 1). It consisted of a pervasive diurnal cycle of seismicity, featuring hundreds of thousands of very small icequakes each day, clustered in time when temperatures were falling. At times of the day when surface temperatures were stable or rising, or at times when there was pervasive standing water on the ice shelf, virtually no seismicity was observed (at frequencies > 1 Hz). The diurnal seismicity cycle was only observed at the wet station and is thus presumed to be caused by processes associated with surface melt/freeze cycles.

The diurnal seismicity cycle had sufficient day-to-day variability to allow its cause to be deduced. For example, we were able to rule out the tidal-flexure process seen elsewhere on ice shelves, most notably at a dry ice-shelf site near the grounding line of Kamb Ice Stream (Hulbe and others, 2016), because days when the seismicity cycle was missing did not show unusually small tidal displacement. (We do not rule out, however, that tidal flexure and other movements do not drive other, smaller-amplitude periods of seismicity on the McMIS. We simply note that tide was not responsible for the conspicuous large diurnal pulsation of interest here.) The physical process that most readily correlates with the seismicity cycle and its variation from day to day was found to be thermally induced bending and fracture of a ‘frozen lid’ that overlies a near-surface partial melt zone.

To test this mechanism, we constructed a simple, idealized thermal model of the ice shelf and simulated the development of subsurface partial melt layers. When forced with AWS-derived weather conditions, the model produced volumes of meltwater, frozen lid thicknesses and thermal bending moments that were consistent with the observed variation in seismicity.

Two types of seismic waves are prominent in our observations. At high frequency, most events are hydroacoustic phases exhibiting a low-frequency cutoff (25–100 Hz). This cutoff is likely induced by the seawater layer between the ice shelf and the seabed which has low sound speed relative to the layers above and below. At low frequency, most events are probably flexural-gravity waves that have a high-frequency cutoff (30 Hz), and that display reverse-dispersion, with high-frequency arriving before low-frequency. These two types of seismic waves follow different population statistics; the high-frequency events feature a Gutenberg--Richter distribution with an intrinsic b-value of 1.2, previously observed in fractured media rupturing at low ambient stress, as is consistent with an ice lid superimposed on a subsurface slush layer. The low-frequency events, whose larger wavelengths would advocate larger source dimensions, expectedly feature a higher b-value, of the order of 3.

We speculate that thermally regulated seismicity may exist on all ice shelves that undergo surface melting and freezing, especially those that have lost their protective, insulating firn layer. We agree, therefore, with previous studies (e.g. Chaput and others, 2018; Podolskiy and Others, 2018) suggesting that seismological methods may be a practical means of monitoring glacier surface conditions over a wide area and through long periods of time. In addition, we note that some seismic signals on ice shelves, such as those associated with surface lake induced hydrofracture, may be masked by seismicity associated with near-surface thermal forcing. The agreement between the numerical model of thermal forcing and seismicity seen in our data suggests it may be possible to use local weather data to separate seismic signals associated with thermal forcing from those due to other phenomena.


We thank the US Antarctic Program, UNAVCO and IRIS/PASSCAL, especially Phillip Chung and Pnina Miller, who have made our study, particularly our presence in the field, possible. We also acknowledge help by Mike Cloutier and the PGC. This work was supported by US National Science Foundation grant PLR-1443126. Seismic data are archived for public distribution at the IRIS Data Management Center (DMC), and is identified by network ID 9G. AWS and GPS data are archived at the US Antarctic Program Data Center (USAP-DC), and are identified by digital object identifier (doi) numbers 10.15784/601106 and 10.15784/601107. AWS data are additionally archived at the University of Wisconsin's Antarctic Meteorological Research Center (AMRC). AFB acknowledges support from a Leverhulme/Newton Trust Early Career Fellowship. GJM acknowledges support from a NASA Earth and Space Science Fellowship (NNX15AN44H). We thank Brad Lipovsky for discussion and valuable suggestions on the analysis presented here. The two anonymous referees and scientific editor who reviewed this manuscript provided extremely good comments and criticism.


DRM, EAO and JL conducted the seismological data analysis. EAO conducted the theoretical analysis presented in the Appendix. AFB, ICW, BG and GJM deployed and recovered the seismographs in the field. All authors were involved in conceiving and writing the manuscript.


Abrahamson, NA and Silva, WJ (1997) Empirical response spectral attenuation relations for shallow crustal earthquakes. Seis. Res. Lett., 68(1), 94127 (doi: 10.1785/gssrl.68.1.94)
Anderson, JG and Lei, Y (1994) Nonparametric description of peak acceleration as a function of magnitude, distance, and site in Guerrero, Mexico. Bull. Seism. Soc. Am., 84, 10031017
Banwell, AF (2017) Ice-shelf stability questioned. Nature, 544(7650), 306307 (doi: 10.1038/544306a)
Banwell, AF and MacAyeal, DR (2015) Ice-shelf fracture due to viscoelastic flexure stress induced by fill/drain cycles of supraglacial lakes. Ant. Sci., 27(6), 587597 (doi: 10.1017/S0954102015000292)
Banwell, AF, MacAyeal, DR and Sergienko, OV (2013) Breakup of the Larsen B Ice Shelf triggered by chain reaction drainage of supraglacial lakes. Geophys. Res. Lett., 40(22), 58725876 (doi: 10.1002/2013GL057694)
Banwell, AF and 5 others (2014) Supraglacial lakes on the Larsen B ice shelf, Antarctica, and at Paakitsoq, West Greenland: a comparative study. Ann. Glaciol., 55(66), 18 (doi: 10.3189/2014AoG66A049)
Banwell, AF and 6 others (2017) Calving and rifting on the McMurdo Ice Shelf, Antarctica. Ann. Glaciol., 58(74), 7887 (doi: 10.1017/aog.2017.12)
Bassis, JN and 7 others (2007) Seismicity and deformation associated with ice-shelf rift propagation. J. Glaciol., 53(183), 523536 (doi: 10.3189/002214307784409207)
Bažant, ZP (1992) Large-scale thermal bending fracture of sea ice plates. J. Geophys. Res. Oceans, 97(C11), 1773917751 (doi: 10.1029/92JC00816)
Bell, RE and 9 others (2017) Antarctic ice shelf potentially stabilized by export of meltwater in surface river. Nature, 544(7650), 344348 (doi: 10.1038/nature22048)
Bromirski, PD and Stephen, RA (2012) Response of the Ross Ice Shelf, Antarctica, to ocean gravity-wave forcing. Ann. Glaciol., 53(60), 163172 (doi: 10.3189/2012AoG60A058)
Bromirski, PD, Sergienko, OV and MacAyeal, DR (2010) Transoceanic infragravity waves impacting Antarctic ice shelves. Geophys. Res. Lett., 37, L02502 (doi: 10.1029/2009GL041488)
Campbell, S, Courville, Z, Sinclair, S and Wilner, J (2017) Brine, englacial structure and basal properties near the terminus of McMurdo Ice Shelf, Antarctica. Ann. Glaciol., 58(74), 111 (doi: 10.1017/aog.2017.26)
Carmichael, JD, Pettit, EC, Hoffman, M, Fountain, A and Hallet, B (2012) Seismic multiplet response triggered by melt at Blood Falls, Taylor Glacier, Antarctica. J. Geophys. Res. Earth Surf., 117, F03004 (doi: 10.1029/2011JF002221)
Cathles, LM IV, Okal, EA and MacAyeal, DR (2009) Seismic observations of sea swell on the floating Ross Ice Shelf, Antarctica. J. Geophys. Res. Earth Surf., 114, F02015 (doi: 10.1029/2007JF000934)
Chaput, J and 10 others (2018) Near-surface environmentally forced changes in the Ross ice Shelf observed with ambient seismic noise. Geophys. Res. Lett., 45 (doi: 10.1029/2018GL079 66 5)
Debenham, F (1965) The genesis of the McMurdo Ice Shelf, Antarctica. J. Glaciol., 5(42), 829832 (doi: 10.3189/S0022143000018888)
Deiz, A and 8 others (2016) Ice shelf structure derived from dispersion curve analysis of ambient seismic noise, Ross Ice Shelf, Antarctica. Geophys. J. Int., 205(2), 785795 (doi: 10.1093/gji/ggw036)
Dudko, YV (1999) Analysis of seismo-acoustic emission from ice fracturing events during SIMI'94, PhD Dissertation, Massachusetts Institute of Technology, Boston, USA, 187 pages
Dunbar, GB, Bertler, NAN and McKay, RM (2009) Sediment flux through the McMurdo Ice Shelf in Windless Bight, Antarctica. Glob. Planet. Change., 69(3), 8793 (doi: 10.1016/j.gloplacha.2009.05.07)
Evans, RJ and Untersteiner, N (1971) Thermal cracks in floating ice sheets. J. Geophys. Res., 76(3), 694703 (doi: 10.1029/JC076i003p00694)
Ewing, M and Crary, AP (1934) Propagation of elastic waves in ice. Part II. Physics, 5(7), 181184 (doi: 10.1063/1.1745249)
Ewing, M and Worzel, JL (1948) Long-range sound transmission. Geol. Soc. Am. Mem., 27, 132 (doi: 10.1130/MEM27-3-p1)
Francis, TJG, Porter, LT and McGrath, JR (1977) Ocean-bottom seismograph observations on the Mid-Atlantic Ridge near latitude 37°N. Geol. Soc. Amer. Bull., 88, 664677
Glasser, NF and Scambos, TA (2008) A structural glaciological analysis of the 2002 Larsen B ice shelf collapse. J. Glaciol., 54(184), 316 (doi: 10.3189/002214308784409017)
Glasser, NF, Goodsell, B, Copland, L and Lawson, W (2006) Debris characteristics and ice-shelf dynamics in the ablation region of the McM Ice Shelf, Antarctica. J. Glaciol., 52(177), 223234 (doi: 10.3189/172756506781828692)
Gutenberg, B and Richter, CF (1954) Seismicity of the Earth and associated phenomena. Princeton University Press, Princeton, N.J., 310 p
Hatherton, T and Evison, FF (1962) A special mechanism for some Antarctic earthquakes. New Zealand J. Geol. Geophys., 5(5), 864873 (doi: 10.1080/00288306.1962.10417642)
Holland, PR and 6 others (2011) The air content of Larsen Ice Shelf. Geophys. Res. Lett., 38, L10503 (doi: 10.1029/2011GL047245)
Hubbard, B and 12 others (2016) Massive subsurface ice formed by refreezing of ice-shelf melt ponds. Nature Comm., 7, 11897 (doi: 10.1038/ncomms11897)
Hulbe, CL and 5 others (2016) Tidal bending and strand cracks at the Kamb Ice Stream grounding line, West Antarctica. J. Glaciol., 62(235), 816824 (doi: 10.1017/jog.2016.74)
Hunkins, K (1960) Seismic studies of sea ice. J. Geophys. Res., 65(10), 34593472 (doi: 10.1029/JZ065i010p03459)
Jarvinen, O and Leppäranta, M (2013) Solar radiation transfer in the surface snow layer in Dronning Maud Land, Antarctica. Polar Science, 7(1), 117 (doi: 10.1016/j.polar.2013.03.002)
Kellogg, DE and Kellogg, TB (1987) Diatoms of the McMurdo Ice Shelf, Antarctica: Implications for sediment and biotic reworking. Paleo. Paleo. Paleo., 60, 7796 (doi: 10.1016/0031-0182(87)90025-3)
Kingslake, J, Ely, JC, Das, I and Bell, RE (2017) Widespread movement of meltwater onto and across Antarctic ice shelves. Nature, 544(7650), 349352 (doi: 10.1038/nature22049)
Klokov, V and Diemand, D (1995) Glaciology of the McMurdo Ice Shelf in the area of air operations. Contributions to Antarctic research IV. Antarc. Res. Ser., 67, 175195 (doi: 10.1002/9781118668207.ch10)
Kuipers Munneke, P and 6 others (2018) Intense winter surface melt on an Antarctic ice shelf. Geophys. Res. Lett., 45, 76157623 (doi: 10.1029/2018GL077899)
Lenaerts, JTM and 12 others (2017) Meltwater produced by wind-albedo interaction stored in an East Antarctic ice shelf. Nature Clim. Chg., 7, 5862 (doi: 10.1038/nclimate3180)
Leppäranta, M, Jarvinen, O and Mattila, O-P (2012) Structure and life cycle of supraglacial lakes in Dronning Maud Land. Ant. Sci., 25(3), 457467 (doi: 10.1017/S0954102012001009)
Lin, J (2017) Seismic Observations on the Melting McMurdo Ice Shelf, Antarctica. MS thesis, University of Chicago, 46 pp
Lipovsky, BP and Dunham, EM (2016) Tremor during ice-stream stick slip. Cryosphere, 10, 385399 (doi: 10.5194/tc-10-385-2016)
MacAyeal, DR (2018) Seismology gets under the skin of the Antarctic Ice Sheet. Geophys. Res. Lett., 45 (doi: 10.1029/2018GL080366)
MacAyeal, DR, Sergienko, OV and Banwell, AF (2015a) A model of viscoelastic ice-shelf flexure. J. Glaciol., 61(228), 635645 (doi: 10.3189/2015JoG14J169)
MacAyeal, DR, Wang, Y and Okal, EA (2015b) Ambient seismic, hydroacoustic, and flexural gravity wave noise on a tabular iceberg. J. Geophys. Res. Earth Surf., 120(2), 200211 (doi: 10.1002/2014JG003250)
Macdonald, G, Banwell, AF and MacAyeal, DR (2018) Seasonal evolution of supraglacial lakes on a floating ice tongue, Petermann Glacier, Greenland. Ann. Glaciol., 59(76pt. 1), 5665 (doi: 10.1017/aog.2018.9)
McNutt, SR (2002) Volcano seismology and monitoring for eruptions. In Lee, WHK, Kanamori, H, Jennings, PC and Kisslinger, C eds. International Handbook of Earthquake and Engineering Seismology, Academic Press, New York, 383406
Milne, AR (1972) Thermal tension cracking in sea ice: a source of under ice noise. J. Geophys. Res., 77(12), 21772192 (doi: 10.1029/JC077i012p02177)
Mogi, K (1962) Magnitude-frequency relation for elastic shocks accompanying fractures of various materials and some related problems in earthquakes (2nd Paper). Bull. Earthq. Res. Inst. Tokyo Univ., 40, 831853
Okal, EA and Kirby, SH (1995) Frequency-moment distribution of deep earthquakes: implications for the seismogenic zone at the bottom of slabs. Phys. Earth Planet. Inter., 92, 169187 (doi: 10.1016/0031-9201(95)03037-8)
Okal, EA and Romanowicz, BA (1994) On the variation of b-value with earthquake size. Phys. Earth Planet. Inter., 87, 5576 (doi: 10.1016/0031-9201(94)90021-3)
Paige, RA (1968) Subsurface melt pools in the McMurdo Ice Shelf, Antarctica. J. Glaciol., 7(51), 511516 (doi: 10.3189/S0022143000020700)
Podolskiy, EA and Walter, F (2016) Cryoseismology. Rev. Geophys., 54(4), 708758 (doi: 10.1002/2016RG000526)
Podolskiy, EA, Fujita, K, Sunako, S, Tsushima, A and Kayastha, RB (2018) Nocturnal thermal fracturing of a Himalayan debris covered glacier revealed by ambient seismic noise. Geophys. Res. Lett., 45, 96999709 (doi: 10.1029/2018GL079653)
Press, F and Ewing, M (1951) Propagation of elastic waves in a floating ice sheet. Trans. of Amer. Geophys. U., 32(5), 673678 (doi: 10.1029/TR032i005p00673)
Rack, W, Haas, C and Langhorne, PJ (2013) Airborne thickness and freeboard measurements over the McMurdo Ice Shelf, Antarctica and implications for ice density. J. Geophys. Res. Oceans, 118(11), 58995907 (doi: 10.1002/2013JC009084)
Rundle, JB (1998) Derivation of the complete Gutenberg-Richter magnitude-frequency relation using the principle of scale invariance. J. Geophys. Res., 94, 1233712342 (doi: 10.1029/93JB02037)
Scambos, T, Hulbe, C, Fahnestock, M and Bohlander, J (2000) The link between climate warming and break-up of ice shelves in the Antarctic Peninsula. J. Glaciol., 46, 516530 (doi: 10.3189/172756500781833043)
Scambos, T and 7 others (2009) Ice shelf disintegration by plate bending and hydro-fracture: satellite observations and model results of the 2008 Wilkins ice shelf break-ups. Earth Plan. Sci. Lett., 280(1–4), 5160 (doi: 10.1016/j.epsl.2008.12.027)
Scholz, CH (1968) The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes. Bull. Seismol. Soc. Amer., 58(1), 399415
Sergienko, OV and MacAyeal, DR (2005) Surface melting on Larsen Ice Shelf, Antarctica. Ann. Glaciol., 40, 215218 (doi: 10.3189/172756405781813474)
Smith, IJ, Langhorn, PJ, Frew, RD, Vennell, R and Haskell, TG (2012) Sea ice growth rates near ice shelves. Cold. Reg. Sci. Technol., 83–84, 5770 (doi: 10.1016/j.coldregions.2012.06.005)
Swithinbank, CWM (1968) Ice movement in the McMurdo Sound area of Antarctica. Int. Symp. Antarctic Glaciological Exploration (ISAGE), Int. Assoc. Sci. Hydrol. SCAR, 86, 472487
Trusel, LD and 6 others (2016) Divergent trajectories of Antarctic surface melt under two twenty-first-century climate scenarios. Nature Geo., 8, 927932 (doi: 10.1038/ngeo2563)
van den Broeke, M (2005) Strong surface melting preceded collapse of Antarctic Peninsula ice shelf. Geophys. Res. Lett., 32, L12815 (doi: 10.1029/2005GL023247)
Walter, F, Deichmann, N and Funk, M (2008) Basal icequakes during changing subglacial water pressures beneath Gornergletscher, Switzerland. J. Glaciol., 54(186), 511521 (doi: 10.3189/002214308785837110)
Winberry, JP, Anandakrishnan, S and Alley, RB (2009) Seismic observations of transient subglacial water-flow beneath MacAyeal Ice Stream, West Antarctica. Geophys. Res. Lett., 36, L11502 (doi: 10.1029/2009GL037730)
Winberry, JP, Anandakrishnan, S, Wiens, DA and Alley, RB (2013) Nucleation and seismic tremor associated with the glacial earthquakes of Whillans Ice Stream, Antarctica. Geophys. Res. Lett., 40, 312315 (doi: 10.1002/grl.50130)
Xie, Y and Farmer, DM (1994) Seismic-acoustic sensing of sea ice wave mechanical properties. J. Geophys. Res., 99(C4), 77717786 (doi: 10.1029/93JC03483)
Yang, TC and Yates, TW (1995) Flexural waves in a floating ice sheet: Modeling and comparison with data. J. Acoust. Soc. of Am., 97(2), 971977 (doi: 1121/1.412076)
Yuan, ATE, McNutt, SR and Harlow, DH (1984) Seismicity and eruptive activity at Fuego Volcano, Guatemala: February 1975 – January 1977. J. Volcan. Geotherm. Res., 21, 277296
Zhan, Z, Tsai, VC, Jackson, JM and Helmberger, D (2014) Ambient noise correlation on the Amery Ice Shelf, East Antarctica. Geophys. J. Int., 196, 17961802 (doi: 10.1093/gji/ggt488)


We consider here the effect of propagation over an unknown (and assumed variable) distance x on the distribution of observed amplitudes, and hence of observed b-values. In very general terms, the decay of a seismic wave with distance has two components, one purely elastic, the other anelastic. The elastic effect is known as geometrical spreading and depends on the nature of the wave: a body wave propagating in a conceptually unlimited, 3-D space will see its amplitude decay as 1/x while a wave propagating in a vertically restricted waveguide (including interface waves such as Rayleigh waves) will feature a $1/\sqrt {x}$ dependence. At very short distances, comparable to the wavelength of the wave, the near field terms of body waves, which scale as 1/x 2, would become prominent. In addition, anelastic attenuation will dampen the amplitude of a seismic wave as

(A1)$$A(x) = A_{{\rm init}}\,e^{-\gamma x}; \quad {\rm with} \quad \gamma = \displaystyle{{\omega} \over {2 Q V}}$$

where ω is the dominant angular frequency of the wave, V its [group] velocity and Q the quality factor of the medium. Because of the exponential nature of Eqn (A1), anelastic attenuation eventually becomes dominant at large distances.

We consider the problem of a population of earthquakes distributed with a parameter b 0, and observed at a variable distance x, the question being ‘Does the population of observed amplitudes satisfy a GR law and, if so, what is the relationship between the observed b 1 at distance x and the original b 0?’

We start with the simple model of a line source where seismicity is distributed along a 1-D axis, with a constant density of source per unit length, μ. The amplitude A observed at some distance x from an event originating with a magnitude m will result from the combination of m and of the attenuation suffered along the path. The combination of geometrical spreading and anelastic attenuation can be represented by a function Φ(x);

(A2)$$\Phi(x) = \Phi_{0} \displaystyle{{\exp {( - \gamma x)}} \over {x}} \quad {\rm or} \quad \Phi(x) = \Phi_0 \displaystyle{{\exp{( - \gamma x)}} \over {\sqrt{x}}}$$

depending on the nature of the wave (body or guided), which controls geometrical spreading. The observed amplitude is then

(A3)$$ A(x,m) = A_0 10^m \cdot \Phi(x) $$

and the apparent ‘magnitude’ M, defined as

(A4)$$ M = \log_{10} \left ( A / A_1 \right ) $$

(where A 0 and A 1 are constants) will simply be

(A5)$$ M = m + \log_{10} \Phi(x) + M_2 $$

with M 2 being yet another constant.

For any given x, the number of events originating between x and x + dx and observed with an apparent magnitude > M will be the number dN of earthquakes with an initial magnitude

(A6)$$ m \gt M - \log_{10} \Phi(x) - M_2 $$


(A7)$$ dN = \mu dx \cdot \displaystyle 10 ^{ C_0 - b_0 \left ( M - \log_{10}\Phi(x) - M_2 \right ) } $$

The total number N of such earthquakes is obtained by simply integrating Eqn (A7) over x. Assuming a simple geometry where the observation point is at the end of the line source of length L,

(A8)$$ N(M) = \mu \int_0^L 10 ^{ C_0 - b_0 \left ( M - \log_{10}\Phi(x) - M_2 \right ) } \cdot {\rm d}x $$

(A9)$$ = \mu 10^{ C_0 - b_0 \left ( M - M_2 \right ) } \cdot \int_0^L \left [ \Phi(x)\right ]^{b_0} \cdot {\rm d} x $$

which takes the form

(A10)$$ N(M) = N_2 \cdot 10^{-b_0 M}; \quad \log_{10}N(M) = C_2 - b_0 M \ \ . $$

It is easily verified that a similar conclusion would be reached if the receiver is outside, or inside, the line source. Its physical interpretation is that, because distance-controlled attenuation is independent of original source size (magnitude m), it affects events of all sizes in the same way and does not perturb the relative distribution of earthquake sizes.

We verified this property by numerically building populations of earthquakes originally satisfying GR laws with various values of b 0, and applying full distance corrections including geometrical spreading. Fig. 23 shows results for a receiver located at the end of the line source (L = 10 km), a 1/x geometrical spreading law, a P-wave velocity in the ice of 3 km s−1, a dominant frequency of 50 Hz and a quality factor of 200; we tested three original b-values, b 0 = 0.5, 1 and 2, respectively. The apparent b-value, b 1, is computed by linearly regressing the cumulative distribution of observed magnitudes. It is clear that the lemma is verified for the small-magnitude portion of the dataset, but that the population of larger events becomes deficient when the magnitude is increased. This is an artifact of the necessity, in our computational exercise, to use an upper limit for the original magnitude m, since we cannot extend it to infinity. When reaching large values of M, the magnitudes m given by Eqn (A5) are undersampled, and we have verified that the domain where this effect becomes significant is a function of the imposed m max, of x, of the parameters defining Φ (e.g. ω and Q), and of the distribution of source size and hence b 0. In retrospect, this apparent limitation in our models is actually reproduced in real life, since the magnitude of the events observed on the ice is expected to be limited upwards by the coherence length over which rupture can develop inside the seismogenic medium.

However, it is important to keep in mind that the development of the elbows in Figs 23 to 25, and the subsequent increases in b 1 are rooted in the presence of a maximum source magnitude m max (which can be physical, as probably is the case for our icequake populations, or artificially imposed by the impossibility of reaching an infinite magnitude in our numerical models). It does not necessarily reflect a genuine change in scaling laws, as for example in the case of the distribution of seismic moments for shallow earthquakes, where an elbow is also observed when the growth of seismic sources is affected by the saturation of fault width at depth (e.g. Okal and Romanowicz, 1994, Fig. 1). In simple terms, an elbow in b 1 does not require one in b 0, if the population of sources is limited to a maximum magnitude m max, but keeps a constant b 0 below that value.

We further tested the robustness of our results with respect to a variation of the geometry of the parameters of the model, with representative examples given on Fig. 24, where we moved the receiver inside the line source (Left), used a much lower level of attenuation (Q = 1000 at f = 10 Hz; Center), or considered a guided wave with a geometrical spreading varying as $1/\sqrt {x}$ (Right), with no significant effect on the main result, i.e., in all cases b 1 ≈ b 0. Finally, in Fig. 25, we explore different geometries, by replacing the line source with a 2-D seismogenic zone of rectangular (Left), square (Center), or circular (Right) shape; once again, we conclude that our results are robust, with b 1 essentially reproducing b 0, especially at low magnitudes.

Fig. 24. Same as Fig. 23 for a receiver located inside the line source (Left); a significantly less attenuating medium (Q = 1000 at f = 10 Hz; Center); and a guided wave with geometrical spreading falling as $1 / \sqrt {x}$ (Right). Note robustness of results.

Fig. 25. Same as Fig. 23 for a receiver located outside a rectangular seismogenic zone (Left); outside a square seismogenic zone (Center); and inside a circular seismogenic zone (Right). Note robustness of results.