A decade of change in the hydraulic connection between an Antarctic epishelf lake and the ocean

Abstract

Observations of the water level in Beaver Lake, an epishelf lake in East Antarctica, show a regular tidal signal that is lagged and attenuated from the tides beneath the adjacent Amery Ice Shelf. The phase lag and amplitude attenuation can be created by a narrow inlet connection between Beaver Lake and the cavity beneath the Amery Ice Shelf. A forced linear damped oscillator is used to determine the inlet dimensions that are required to produce the observed phase lag and amplitude attenuation. The model shows that the observations are consistent with a tidal flow that is restricted by the drag created by flow in the narrow inlet. Analysis shows that the phase lag and amplitude attenuation of the tides in Beaver Lake has increased over the years 1991-2002, probably due to a thickening of the overlying ice shelf. The response is sensitive to subtle variations in the dimensions of the inlet.

1. Introduction

Epishelf lakes are freshwater, tidal lakes that are connected to the open ocean via paths beneath floating ice shelves. Epishelf lakes are located between rocky ice-free areas and floating ice shelves in the polar regions and have been found in Arctic Canada (Vincent and others, 2001) and in areas around the Antarctic continent (Gibson and Andersen,. The existence of epishelf lakes is dependent on freshwater resupply and the continued embayment by a floating ice shelf. The unique physical position of epishelf lakes makes them especially sensitive to changes in climate (e.g. Vincent and others, 2001; Veillette and others, 2008). The collapse of the embaying ice shelf can have catastrophic consequences for the epishelf lake − the abrupt collapse of the Ward Hunt Ice Shelf-Disraeli Fjord, leading to the loss of the adjoining epishelf lake, has been well documented (Vincent and others, 2001; Mueller and others,. Epishelf lakes can host unusual communities of marine and freshwater biota (e.g. Laybourn-Parry and others, 2001) and are unique places for the deposition of plankton species useful for palaeoclimate records (e.g. Smith and others, 2006).

One of the largest known Antarctic epishelf lakes is Beaver Lake (BL), located at 70°48’S, 68°15'E in the Amery Ice Shelf (AIS) region of East Antarctica (Fig. 1). BL is classified as a type 1 epishelf lake; the freshwater layer floats on saline marine water, and the depth of the freshwater layer (~200 m; Gibson and Andersen, 2002) is controlled by the thickness of the ice shelf. There have been many studies to try to assess the physical (e.g. Bardin and others, 1990; Gibson and Andersen, 2002) and biological (e.g. Laybourn-Parry and others, 2001) conditions of the lake waters.

Fig. 1. The location of Beaver Lake and the Amery Ice Shelf in East Antarctica (black box in inset), adapted from Galton-Fenzi and others (2008).

The connection with the ocean makes the BL tide record an important input to AIS research. The rocky shoreline provides a stable environment for vertically referencing sea level. However, studies have shown that the BL tides lag the tides beneath the AIS by ~1 hour (e.g. King and others, 2000; Fricker and others, 2002). Hydrodynamic models of the AIS cavity, which include BL, cannot account for the phase lag and amplitude attenuation of the tides and suggest unresolved local topographical features as the explanation (Padman and others, 2002; Hemer and others, 2006; Maraldi and others, 2007).

Here we show that the linear response of an almost isolated lake basin to a tide in the ocean beneath an ice shelf can be used to describe the phase lag and amplitude attenuation. Furthermore, we use three separate sets of tidal observations, made over an 11 year period, which show that the phase lag and amplitude attenuation has increased. The physical response of the inlet that is required to produce the observed phase lag and amplitude attenuation is addressed analytically and numerically. The inlet dimensions can be inferred from some basic estimates of the lake area and tidal observations. The surface area within the lake is used to constrain the solution and is estimated from synthetic aperture radar (SAR) interferograms. Analysis using constraints from observations suggests that the inlet has narrowed due to a decrease in the net basal melt at the base of the ice shelf overlying the inlet. Note that the lake referred to in this study encompasses the entire area within the solid line in Figure 2 that is south of ~70°40′S. The patches of open water can be seen in Figure 2a as smooth areas, and the textured region is the overlying ice shelf.

Fig. 2. (a) Beaver Lake, showing the locations of the tide gauges used in the 1990/91 and 1997/98 seasons as a filled circle (70º48′ S, 68º10′ E) and the GPS receiver used in 2002 as a cross (70º47′ S, 68º13′ E) overlaid on a Landsat image from 2005. (b) SAR interferogram showing an outline of Beaver Lake (solid line) and the location of the grounded region confining the inlet (dashed line), generated from two pairs of European Remote-sensing Satellite ERS-1 and ERS-2 SAR data from the tandem mission acquired in February and March 1996.

2. Theory

Tidal elevations in a semi-enclosed basin (such as BL), connected to the sea through a narrow inlet, may be approximated by a linear damped oscillator. It may be assumed that the sea level is nearly uniform within the basin and has a constant slope within the inlet (Fig. 3). The response in the lake is determined by the near-resonance strength of the tidal forcing and the strength of the damping. Here damping is controlled by two main processes: (1) drag due to frictional stress on the inlet walls and (2) drag due to flow separation, created by the difference between the dynamic pressure and the total pressure difference driving the flow through the inlet. These two processes lead to the same linear quadratic damping law.

Fig. 3. Schematic of a lake connected to the ocean via an inlet. A is the horizontal surface area of the basin, η _b is the surface elevation within the basin (assumed spatially constant), η _b is the surface elevation in the open ocean, u _c is the horizontal velocity in the channel (assumed spatially constant) and L, h and W are the length, depth and width of the channel, respectively (assumed constant).

Conservation of volume for the basin yields

(1)

where A is the horizontal surface area of the basin, η _b is the surface elevation within the basin (assumed spatially constant), t is time, u _c is the horizontal velocity in the inlet (assumed spatially constant) and h and W are the depth and width of the inlet, respectively (assumed constant).

Any change in surface elevation in the basin is due to extra volume entering through the inlet. Conservation of momentum for the inlet, without rotation or inertia, yields

(2)

where g is acceleration due to gravity, η _c is the (longitudinally varying) surface elevation in the inlet, τ is the frictional stress on the walls of the inlet, P is a dynamic pressure term parameterizing flow separation at the ends of the inlet and ρ is the water density. L is the length of the inlet and R is the hydraulic radius, which is the ratio of the cross-sectional area of the inlet to its wetted perimeter (noting that the inlet is enclosed on the upper surface by an ice shelf):

(3)

The flow is driven by the pressure difference, ρg(η _b – η _°) where η _º is the surface elevation in the open ocean. The friction term, which may be linearized, is proportional to the typical current magnitude, yielding τ/ρ = C _f u _b u _c, where C _f ≈ 0.0025 is an empirical quadratic drag coefficient due to skin friction on the walls of the inlet and u _b is a constant depth-averaged velocity of approximately the same magnitude as u _c. Similarly, the dynamic pressure term may be linearized, yielding P = _ρ C _p u _b u _c, where C _p is the coefficient of head loss due to flow separation. Assuming a finite-difference approximation, ∂η _c/∂x ≈ (η _° — η_b)/L, yields

(4)

The value of C _p has been found empirically to range from 0.97 for idealized laboratory experiments (Terra and others, 2005) to 1.5 for geophysical settings (summarized by Miles and Lee, 1975). Here we assume C _p = 1.

Combining Eqns (1) and (4) and incorporating the linear damping terms yields:

(5)

Where

(6a)

(6b)

For harmonic solutions of the form η_b ∞ e^(iωt) and η_° ∞ e^(iωt)

(7)

where ω is the angular frequency of the tidal forcing, which has an amplitude ratio, |η_b/η_°|, and phase lag, ϕ, given by

(8a)

(8b)

where ϕ is positive when η_b lags η_°. Note that this system has a single resonance in the region of the fundamental or Helmholtz mode. The sharpness of the Helmholtz mode depends on the quality of the resonance, which is summarized by the value of the ‘Q-factor’, given by (a high Q indicates a sharp resonance).

By determining |η_b/η_°| and ϕ from tidal analysis (which is done for BL in Section 3.2) then we may solve Eqns (8a) and (8b) for α and β:

(9a)

(9b)

Approximate solutions can be found by assuming u _b ≈ |u _c| which, from Eqn (1), is:

(10)

Estimating any one of the three unknown inlet parameters (h, L or W) then the remaining two can be found using Eqns (6a), (6b) and (10). This solution also requires estimates of a and β, using Eqns (9a) and (9b), from tidal observations and an estimate of the surface area of the lake, A.

3. Measurements

3.1. Geometry

Defining the dimensions of BL, particularly the inlet geometry, is problematic as there are no in situ observations. For estimates of the size of the lake we rely on SAR observations. The SAR interferogram (Fig. 2b) is a doubledifference product generated from two pairs of ERS-1 and ERS-2 SAR data from the tandem mission acquired in February and March 1996. SAR interferometry is a well-established technique (e.g. Joughin and others? 1999; Young and Hyland, 2002). The interferogram fringe pattern represents contributions from a combination of surface topography and vertical motion of the ice at the zone located between grounded ice and floating ice. Field observations suggest that much of the lake is covered by the ice shelf, the rest being covered for most of the year by fast ice. Tide cracks that have been seen around the lake edge suggest the ice shelf within BL is mostly free-floating. The hinge zone around the southern end of the lake (derived from SAR) resembles the lake boundary shown in Figure 2a. Most of the ice-covered region in BL is free-floating (Fig. 2b). The free surface area of the lake can therefore be estimated directly from maps of the area. However, there is some uncertainty in discriminating between features near the inlet. The area contained between the dashed and solid lines in Figure 2b can be interpreted in two ways: (1) as permanently grounded and forming a long and narrow inlet (solid line) or (2) as ephemerally grounded and forming a wide and short inlet (dashed line). The effective lake area, A, is ~800km² and the inlet width is 3−8km. These measurements are used in the analysis to determine the effective inlet dimensions that predict the observations.

3.2. Tides

Historic observations of the BL water level during 1958 revealed a regular tidal signal. The observations, though too short for our purposes, showed that a connection between BL and the ocean beneath the AIS did exist (McCleod, 1959). Further tidal observations were made by the Australian National Tidal Centre using a bottom-mounted pressure gauge over two periods: 26 December 1990 to 18 January 1991 (22.85 days; N.O. Ward, unpublished information) and 1 December 1997 to 1 March 1998 (90 days; J. Hyslop, unpublished information). The gauges were deployed within 500 m of each other, at an approximate location of 70°48' S, 68°10' E (Fig. 2a, solid circle). GPS data were also collected at a station located on the fast ice (70°47′ S, 68° 13'E (Fig. 2a, +) from 14 January to 9 February 2002 (26 days). The GPS position was computed every 5 min.

Tidal constant analysis (amplitude and phase) was done using the T_TIDE harmonic analysis software (Pawlowicz and others, 2 002). The analysis concentrated on the major diurnal (K₁ and O₁) and semi-diurnal (S₂ and M₂) constituents. Due to the shortness of the BL records, the analysis separated P₁ from K₁ and K₂ from S₂ using the sea-level record from Mawson (67°36′ S, 62°52′ E), the nearest long-term tide observations to BL. This analysis uses the well-known response method of Munk and Cartwright (1966) to separate the tidal constituents. The tidal observations are presented in Table 1.

Table 1. Observed tides in Beaver Lake. The amplitude, h (m), with error, ∊ _h , and phase, g, with error, ∊ _g , for the major diurnal tides K₁ and O₁ and the major semi-diurnal tides M₂ and S₂. These four constituents contribute ~80% of the tide in Beaver Lake

Observations of the tides beneath the AIS are not available for the same period as the observations in BL. Instead, the Circum-Antarctic Tidal Simulation (CATS02.01) developed by Padman and others (2002) was used to predict the major diurnal (O₁, K₁) and semi-diurnal (M₂, S₂) tides. CATS02.01 has been assessed in three separate studies of the barotropic tides of the AIS region (King and others, 2005; Hemer and others, 2006; Maraldi and others, 2007). The CATS02.01 simulations were found to have good agreement with observations in the open ocean, though tend to overpredict the amplitude of the S₂ constituent. The tides were predicted for the closest model cell to BL (70° S, 70° E). The tidal predictions are presented in Table 2.

Table 2. Simulated tides in the ocean beneath the Amery Ice Shelf. The amplitude, h (m), with error, ∊_h , and phase, g, with error, ∊_g , for the major diurnal tides K₁ and O₁ and the major semi-diurnal tides M₂ and S₂

Observations show that the major tidal semi-diurnal (M₂ and S₂) and diurnal (O₁ and K₁) constituents explain up to 80% of the tidal variance. The amplitude ratio and phase differences between the observed tides in BL and the simulated tides beneath the AIS are presented in Table 3. The phase lag and attenuated amplitude exists for both the semi-diurnal and diurnal tidal constituents, with the most obvious difference between years in the phase difference of the diurnal tidal constituents. An amplitude attenuation is seen for the remaining constituents. All constituents but the K₁ show an increasing phase lag and amplitude attenuation with time. The decreased amplitude attenuation of K₁ is not statistically significant.

Table 3. The amplitude ratio, |η_b/η_°|, and phase difference, ϕ, and associated errors, and ∊_|ηb/η°| and ∊_ϕ , between observations within BL and simulated tides beneath the AIS. Major diurnal tides K₁ and O₁ and major semi-diurnal tides M₂ and S₂

The phase lags are largest for the diurnal constituents and show a statistically significant increase over time. However, the analysis shows that the S₂ constituent is amplified in BL, which we attribute to poor representation in the tide model. For the analysis of the analytic model we utilize the phase lag and amplitude attenuation of the M₂ tidal constituent, as it shows the statistically significant changes and is the most robustly predicted.

4. Discussion

BL can be approximated by the linear model described above, where the source of damping is due to flow separation and drag on the inlet perimeter. This then allows the inlet cross-sectional area and inlet length to be estimated independently from tidal observations. As an example, the M₂ tidal constituent is used here to calculate the solutions, as it shows a significant phase lag and amplitude attenuation in BL and is accurately predicted in the AIS cavity. Taking the appropriate phase lag and amplitude attenuation for M₂ from Table 3, we can find appropriate values of h and L for a known W by solving Eqns (6a), (6b) and (10), as

(11a)

(11b)

Where

Solutions using these equations, for a range of W, are shown in Figure 4. We can see that realistic values for the inlet dimensions exist for all years. The height, particularly, which is the least-known parameter, is predicted within a reasonable range. The SAR interferogram shows that much of the inlet is freely floating. Thus, the clearance between the upper ice surface and sea floor should be larger than the range of sea-level movement over a full tidal cycle, which can be ~1-2 m. Interpreting the SAR interferogram as showing a region of grounded ice (Fig. 2b), then by using characteristic parameters from 1997/98, the appropriate inlet length can be calculated as ~650m for a mean inlet height of ~10.2 m.

Fig. 4. The scaling of the inlet height, h, and length, L, with the width, W, using the characteristic parameters shown in Table 4 for 1990/91 (solid line), 1997/98 (dashed line) and 2002 (dashed-dotted line).

The SAR interferogram suggests that the inlet length should be much larger than this, at ~20km. The likely reason for the difference is that the constriction, which acts to attenuate the tides, is only important for a short length, possibly due to the sea floor becoming shallower in the region of the inlet. Note that, without knowledge of the sea floor and ice draft, the short constriction could either be created by a shoaling seabed or by icy protuberances that could be hanging from the underside of the ice shelf. We consider the latter to be unlikely, as there is no evidence for such large changes in the surface elevation.

The model shows that the increased phase lag and amplitude attenuation over the 11 year period is caused by an inlet that is becoming both narrower and shorter. The effect of the constricting inlet can be seen in Q (Table 4). The Q-factor shows that the changing dimensions of the inlet can be reflected in a deterioration in the quality of the lowest resonant period. The model also shows that the actual effective inlet length is relatively short compared with the length of BL. That is, L is only important over the region of maximum constriction (see Fig. 5).

Table 4. Characteristic parameters α and β determined from Eqns (9a) and (9b), using M₂ tidal observations. From these and A = 800 km², ω = 1.4053 × 10^–4s–¹ and η_b, Helmholtz period 2π/ω_H and Q-factor Q are obtained. The inlet width, W, height, h, and length, L, are shown in Figure 4

Fig. 5. Schematic of an epishelf lake and an overlying ice shelf that creates a short (small L) inlet. The fresh water can be supplied from a combination of surface runoff and basal melt from the embaying ice shelf. The extent of the lake is shown by the horizontal line.

The condition required to create a narrower inlet can arise from a change in the mass balance of Charybdis Glacier, which flows into BL. This can happen in two ways: (1) the flow rate of Charybdis Glacier into BL can increase or (2) the melting at the base of Charybdis Glacier that overlies the inlet can decrease.

A comparison of Landsat data from 1973 with data from 2005 shows no obvious change in the southern extent of the ice in the BL basin. Tracking features in the ice shelf from 27 December 1999 to 17 November 2002 (at 70°48'S, 68°16'E) show a southward movement (195- 205°True North (TN)) of the ice of ~62±8m a~¹. This agrees with GPS measurements which show movement of the fast ice in a direction of ~200° TN at ~82 m a^-1.

Thus, the overall extent of Charybdis Glacier within BL has not changed over the period of sea-level observations. The fact that the Landsat observations show no movement of the overlying ice shelf into BL suggests that the inlet may be becoming more constricted due to changes in the basal conditions. This may be due to a combination of refreezing and less melting at the ice-shelf base, leading to a thickening of the overlying ice shelf. Two mechanisms that can cause local cooling of the waters in contact with the base of the ice shelf are: (1) a decrease in the relatively warm freshwater runoff (compared to the freezing point of sea water) into BL and/or (2) a cooling of the ocean water that comes from beneath the AIS.

The thickening discussed here is broadly consistent with observations of an increase in the surface elevation of the AIS from 1996 to 2003 (King and others, 2009). The increase in surface elevation is interpreted as suggesting that the AIS thickened from 1992 to 2003 but otherwise has been relatively stable on multi-decadal timescales. A possible consequence of the thickening is that BL might eventually become isolated from the ocean beneath the AIS. However, King and others (2009) suggest a more recent thinning trend, over 2003-07. Our observations do not span this epoch but we can suppose that if the activity of the ice shelf overlying the inlet is consistent with the observations on the AIS, we can expect that the tidal lag and attenuation has decreased. Both of these scenarios have implications for the communities of marine and freshwater biota that live within the lake.

5. Conclusions

A phase lag and amplitude attenuation of sea surface has been observed in the tides between BL and the ocean beneath the AIS. We have shown that this phase lag can be recreated using a simple model of a linear damped oscillator. We have analysed three sets of records spanning 11 years and observe that the lag has increased with time. An increased phase lag between BL and the ocean beneath the AIS from 1991 to 2002 can be described by an increase in the inlet friction which suggests both that the inlet has become more constricted and the effective inlet length has decreased.

The thickening is consistent with observations of thickening on the AIS (King and others, 2009). Observations of the water level inside BL and in the adjacent AIS, when analysed with the model outlined here, may provide a simple method to monitor long-term changes to the ice shelf. The results suggest that the region of grounded ice in the middle of the inlet has an important effect on the hydraulic connection with the ocean beneath the AIS, and this has been included in recent work describing the cavity geometry (Galton-Fenzi and others, 2008). However, the solutions are sensitive to small changes in the inlet dimensions. Hydrodynamic models in this region will need a high resolution (of the order of metres) to correctly resolve the tidal response inside BL.