## References

Albrecht, T and Levermann, A (2012) Fracture field for large-scale ice dynamics. Journal of Glaciology 58(207), 165–176. doi: 10.3189/2012JoG11J191

Arbic, BK, Mitrovica, JX, MacAyeal, DR and Milne, GA (2008) On the factors behind large Labrador Sea tides during the last glacial cycle and the potential implications for Heinrich events. Paleoceanography 23, PA3211. doi: 10.1029/2007PA001573

Arbic, BK, Karsten, RH and Garrett, C (2009) On tidal resonance in the global ocean and the back-effect of coastal tides upon open-ocean tides. Atmosphere-Ocean 47(4), 239–266. doi: 10.3137/OC311.2009

Aubrey, DG and Speer, PE (1985) A study of non-linear tidal propagation in shallow inlet/estuarine systems Part I: observations. Esturarine, Coastal and Shelf Science 21(2), 185–205. doi: 10.1016/0272-7714(85)90096-4

Banwell, AF and MacAyeal, DR (2015) Ice-shelf fracture due to viscoelastic flexure stress induced by fill/drain cycles of supraglacial lakes. Antarctic Science 27(6), 587–597. 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. Geophysical Research Letters 40(22), 5872–5876. doi: 10.1002/2013GL057694

Banwell, AF, Willis, IC, Macdonald, GJ, Goodsell, B and MacAyeal, DR (2017) Calving and rifting on the McMurdo Ice Shelf, Antarctica. Annals of Glaciology 58(75pt1), 78–87. doi: 10.1017/aog.2017.12

Banwell, AF, Willis, IC, Macdonald, GJ, Goodsell, B and MacAyeal, DR (2019) Direct measurements of ice-shelf flexure caused by surface meltwater ponding and drainage. Nature Communications 10(1), 730. doi: 10.1038/s41467-019-08522-5

Banwell, AFand 5 others (2014) Supraglacial lakes on the Larsen B ice shelf, Antarctica, and at Paakitsoq, West Greenland: a comparative study. Annals of Glaciology 55(66), 1–8. doi: 10.3189/2014AoG66A049

Brunt, KM and MacAyeal, DR (2014) Tidal modulation of ice-shelf flow: a viscous model of the Ross Ice Shelf. Journal of Glaciology 60(221), 500–508. doi: 10.3189/2014JoG13J203

Fricker, HA and Padman, L (2006) Ice-shelf grounding-zone structure from ICESat laser altimetry. Geophysical Research Letters 33, L15502. doi: 10.1029/2006GL026907

Georgas, N (2012) Large seasonal modulation of tides due to ice cover friction in a midlatitude estuary. Journal of Physical Oceanography 42(3), 352–369. doi: 10.1175/JPO-D-11-063.1

Glasser, N, Goodsell, B, Copland, L and Lawson, W (2006) Debris characteristics and ice-shelf dynamics in the ablation region of the McMurdo Ice Shelf, Antarctica. Journal of Glaciology 52(177), 223–234. doi: 10.3189/172756506781828692

Glasser, NF and Scambos, TA (2008) A structural glaciological analysis of the 2002 Larsen B ice-shelf collapse. Journal of Glaciology 54(184), 3–16. doi: 10.3189/002214308784409017

Gleason, CJand 8 others (2016) Characterizing supraglacial meltwater channel hydraulics on the Greenland Ice Sheet from in situ observations. Earth Surface Processes and Landforms 41, 2111–2122. doi: 10.1002/esp.3977

Henderschott, MC (1980) Long waves and ocean tides. In Warren BA and Wunsch C (eds), *Evolution of Physical Oceanography. Scientific Surveys in Honor of Henry Stommel*. Cambridge: MIT Press, pp. 292–341.

Hutter, K, Wang, Y and Chubarenko, IP (2014) Physics of Lakes. Berlin: Springer International. 1685 pp.

King, MAand 6 others (2011) Ocean tides in the Weddell Sea: new observations on the Filchner-Ronne and Larsen C ice Shelves and model validation. Journal of Geophysical Research Oceans 116, C06006. doi: 10.1029/2011JC006949

Kingslake, J, Ely, JC, Das, I and Bell, RE (2017) Widespread movement of meltwater onto and across Antarctic ice shelves. Nature 544(7650), 349–352. doi: 10.1038/nature22049

Kingslake, J, Ng, F and Sole, A (2015) Modelling channelized surface drainage of supraglacial lakes. Journal of Glaciology 61(225), 185–199. doi: 10.3189/2015JoG14J158

Knudsen, P, Bingham, R, Andersen, O and Rio, M-H (2011) A global mean dynamic topography and ocean circulation estimation using a preliminary GOCE gravity model. Journal of Geodesy 85(11), 861–879. doi: 10.1007/s00190-011-0485-8

Lenaerts, JTMand 12 others (2017) Meltwater produced by wind-albedo interaction stored in an East Antarctic ice shelf. Nature Geoscience 7(1), 58–62. doi: 10.1038/nclimate3180

MacAyeal, DR (1984) Numerical simulations of the Ross Sea tides. Journal of Geophysical Research Oceans 89(C1), 607–615. doi: 10.1029/JC089iC01p00607

MacAyeal, DR and Sergienko, OV (2013) The flexural dynamics of melting ice shelves. Annals of Glaciology 54(63), 1–10. doi: 10.3189/2013AoG63A256

MacAyeal, DR, Sergienko, OV and Banwell, AF (2015) A model of viscoelastic ice-shelf flexure. Journal of Glaciology 61(228), 635–645. doi: 10.3189/2015JoG14J169

Macdonald, GJand 5 others (2019) Formation of pedestalled, relict lakes on the McMurdo Ice Shelf, Antarctica. Journal of Glaciology 65(250), 337–343. doi: 10.1017/jog.2019.17

Mortimer, CH and Fee, EJ (1976) Free surface oscillations and tides of lakes Michigan and Superior. Philosophical Transactions of the Royal Society of London, Series A 281(1299), 1–61. doi: 10.1098/rsta.1976.0020

Mueller, RDand 5 others (2012) Impact of tide-topography interactions on basal melting of Larsen C Ice Shelf, Antarctica. Journal of Geophysical Research Oceans 117, C05005. doi: 10.1029/2011JC007263

Padman, L, Erofeeva, S and Joughin, I (2003b) Tides of the Ross Sea and Ross Ice Shelf cavity. Antarctic Science 15(1), 31–40. doi: 10.1017/S0954102003001032

Padman, L, Fricker, HA, Coleman, R, Howard, S and Erofeeva, L (2003a) A new tide model for Antarctic ice shelves and seas. Annals of Glaciology 34, 247–254. doi: 10.3189/172756402781817752

Padman, L, Siefgried, MR and Fricker, HA (2018) Ocean tide influences on the Antarctic and Greenland ice sheets. Reviews of Geophysics 56, 142–184. doi: 10.1002/2016RG000546

Pedlosky, J (1987) Geophysical Fluid Dynamics, 2nd Edn.New York: Springer-Verlag. 710 pp. doi: 10.1007/978-1-4612-4650-3

Platzman, GW (1978) Normal modes of the world ocean. Part I. Design of a finite-element barotropic model. Journal of Physical Oceanography 8(3), 323–343. doi: 10.1175/1520-0485(1978)008<0323:NMOTWO>2.0.CO;2

Platzman, GW, Curtis, GA, Hansen, KA and Slater, RD (1981) Normal modes of the World Ocean. Part II: description of modes in the period range 8 to 80 hours. Journal of Physical Oceanography 11(5), 579–603. doi: 10.1175/1520-0485(1981)011<0579:NMOTWO>2.0.CO;2

Rack, W, Haas, C and Langhorne, PJ (2013) Airborne thickness and freeboard measurements over the McMurdo Ice Shelf, Antarctica, and implications for ice density. Journal of Geophysical Research Oceans 118(11), 5899–5907. doi: 10.1002/2013JC009084

Robel, A and Banwell, AF (2019) A speed limit on ice shelf collapse through hydrofracture. Geophysical Research Letters 46(21), 1292–12100. doi: 10.1029/2019GL084397

Robertson, R, Padman, L and Egbert, GD (1985) Tides in the Weddell Sea. In Jacobs SS and Weiss RF (eds), *Ocean, Ice, and Atmosphere, Interactions at the Antarctic Continental Margin*. Antarctic Research Series, Vol. 75. American Geophysical Union: Washington, DC, pp. 341–369. doi: 10.1029/AR075p0341

Rosier, Sand 5 others (2017) On the interpretation of ice-shelf flexure measurements. Journal of Glaciology 63(241), 783–791. doi: 10.1017/jog.2017.44

Rosier, SHR, Green, JAM, Scourse, JD and Winkelmann, R (2014) Modeling Antarctic tides in response to ice shelf thinning and retreat. Journal of Geophysical Research Oceans 119(1), 87–97. doi: 10.1002/2013JC009240

Scambos, T, Hulbe, C and Fahnestock, M (2003) Climate-induced ice shelf disintegration in the Antarctic Peninsula. In Domack EW, Burnett A, Leventer A, Conley B, Kirby M and Bindschadler R (eds), *Antarctic Peninsula Climate Variability: A Historical and Paleoenvironmental Perspective*. Washington, DC: American Geophysical Union, pp. 79–92. doi: 10.1029/AR079p0079.

Smith, LCand 22 others (2017) Direct measurements of meltwater runoff on the Greenland Ice Sheet surface. Proceedings of the National Academy of Sciences 114(50), E10622–E10631. doi: 10.1073/pnas.1707743114

Vaughan, DG (1995) Tidal flexure at ice shelf margins. Journal of Geophysical Research Solid Earth 100(B4), 6213–6224. doi: 10.1029/94JB02467

Walker, RTand 5 others (2013) Ice-shelf tidal flexure and subglacial pressure variations. Earth and Planetary Science Letters 361, 422–428. doi: 10.1016/j.epsl.2012.11.008

Wild, CT, Marsh, OJ and Rack, W (2018) Unraveling InSAR observed Antarctic ice-shelf flexure using 2-D elastic and viscoelastic modeling. Frontiers of Earth Science 6, 1–11. doi: 10.3389/feart.2018.00028

Zaron, ED (2018) Ocean and ice shelf tides from CryoSat-2 altimetry. Journal of Physical Oceanography 48(4), 975–993. doi: 10.1175/JPO-D-17-0247.1

## APPENDIX A. Derivation of model equations

The shallow-water equations applicable to the surface meltwater layer are written

where *η*, *ζ* and *T* are the water-depth perturbation, the elastic ice-shelf deflexion and the ocean tide, respectively, where

and where ${\bar \cdot}$ denotes the long-term average operator, *u* and *v* are the *x* and *y* components of meltwater velocity ${\bf u}$ (assumed independent of *z*, the vertical coordinate), *g* is the acceleration of gravity, and *M* is the Manning roughness coefficient taken to range from 0.01 to 0.15 s m^{1/3} (Smith and others, 2017). In writing the above equations, we have chosen to disregard Coriolis forces to simplify the analysis here. This disregard is justified by the fact that the horizontal span of meltwater regimes on ice shelves is order 1–10 km in scale, and rotation effects are commonly disregarded in such small water bodies (Mortimer and Fee, 1976).

Our first simplification is to linearize the equations by disregarding the ${\bf u} \cdot \nabla$ terms, replacing *h* with ${\bar h}$, and writing the friction in terms of a frictional timescale *τ*:

Taking *h* ~ 1 m and $\left \vert {\bf u} \right \vert \sim 0.01$ m s^{−1}, we have

Equations (A1) and (A2) are now written,

We now operate on the first of the above equations with the operator

and use the second equation to replace the result of the operator acting on ${\bf q}$ in the right-hand side of the resulting equation to get:

We now non-dimensionalize the variables using scales *L* ~ 10^{3} m for *x* and *y*, *h* _{o} ~ 1 m for ${\bar h}$, *T* _{o} ~ 1 m for *η*, *ζ* and *T*, and $\sigma = \lpar {2 \pi }/{24}\rpar \, {\rm h}^{-1}$ for ∂/∂*t* terms. The result, after a bit of algebra, is

where the variables are now scaled in dimensionless form. We now evaluate the two parameters appearing on the left-hand side of the above equation:

We thus disregard the acceleration term involving ∂^{2}*η*/∂*t* ^{2} term on the left-hand side of the above equation while retaining the ∂*η*/∂*t* term, giving as the final simplified result:

where we have dropped the ${\bar \cdot}$ notation over *h*, and which is a diffusion-type equation for the variable *η*. This constitutes the derivation of Eqn (A14).

The equation governing ice-shelf flexure ζ(*x*, *y*, *t*) treats the ice shelf as a thin elastic plate lying above an infinite half-space ocean (MacAyeal and others, 2015):

where *ρ* _{sw} = 1028 kg m^{−3} is the density of seawater, *ρ* _{w} = 1000 kg m^{−3} is the density of meltwater on top of the ice shelf, *M* _{xx}, *M* _{xy} and *M* _{yy} are the bending moments which are related to the curvature of the flexure by the following expressions:

The variable *H* is the effective ice-shelf thickness (with regard to flexure) and is taken to be constant at 30 m which is consistent with what is observed near Rift Tip Lake (Rack and others, 2013), *E* = 5 GPa is the Young's modulus for ice, and is taken to be in the middle of the observed range of variation for natural ice (1–10 GPa), and *μ* = 1/3 is the Poisson ratio for ice. Elastic parameters and ice-shelf thickness influence the resultant flexure in a manner that is captured by a single parameter $\cal F$, the elastic rigidity:

which shows that the resistance to flexure increases with the cube of the ice thickness. This parameter is also involved in determining the horizontal length scale of ice-shelf flexure effects, and hence what size the surface meltwater lake must be to be significantly influenced by flexure, and *vice versa*.

## APPENDIX B. Sensitivity of supraglacial meltwater response tidal tilt

Although our study concludes that tidal tilt was not the likely cause of diurnal meltwater-depth variations on the McMIS, our study did not rule out the possibility that tidal tilt could be a significant factor in driving surface meltwater oscillations in some circumstances (i.e. surface meltwater regimes with large horizontal scale in regions of high tidal amplitude variation). Even though tidal-tilt-driven meltwater phenomena may ultimately prove to be negligible, it is important that future investigation of its possible workings are appraised so as to avoid negative confirmation bias. We thus examine the response to tidal tilt in a more comprehensive manner in this appendix.

Equations (A14)–(A17) used in this study represent a highly simplified (linearized, simplified friction and boundary conditions) model of surface meltwater behavior in the ice-shelf/ocean system, and depends on three non-dimensional parameters:

where Δ*ρ* = *ρ* _{sw} − *ρ* _{w} ~ 28 kg m^{−3} is the density difference between salty seawater and meltwater. The model equations are now examined in light of the possible variation of the above parameters to determine representative behaviors that may exist under the great variety of conditions expected on an ice shelf.

### Regime I. Friction-free, rigid ice shelf

If τ ≫ 10^{3} s and *L* is appropriately small (<10 km), we expect *κ* > 10^{2}, and if we also consider the ice-shelf to be sufficiently stiff, with ${\cal F} \gg 1$ (e.g. by requiring *H* to be 100 m or greater), the governing equations (in dimensional form) reduce to

The first of the above equations is twice integrable to yield,

where $\cal A$ is the region of each closed, internally connected meltwater system. The integral in the above solution is the constant of integration (the linear term depending on *x* and *y* coming from the first integral of Eqn B10 must be zero to satisfy zero gradient, and thus zero meltwater flux, at the boundaries of $\cal A$ denoted by $\delta \cal A$) that is required to ensure that

which is a statement of mass conservation within each closed meltwater system. Under circumstances where the meltwater system may be only partially closed, e.g. when a spillway exists for certain ranges of water level in certain places, the above constraint may not hold; but we shall not consider this while investigating tidal-tilt phenomena.

The solution given by Eqn (B12) shows that the meltwater movement is such as to ensure that the surface of the meltwater layer is flat through all phases of the tidal cycle. The very large value of the effective diffusivity ensures that meltwater transport is able to keep up with the gentle tilts of the ocean tide and supply extra water to those parts of the meltwater system that would otherwise have a different surface level.

Under circumstances where we have a local linearization of the tidal range, e.g. where, as an illustration,

we would then have,

Under most circumstances, where *L* is on the order of kilometers at most, the magnitude of *T*′ is probably at most in the millimeter or centimeter range. This means that when the ice shelf is stiff, but the meltwater layer highly mobile, the sympathetic tide driven by *T* is of minimal amplitude. A schematic diagram of this low-friction, low-ice-shelf rigidity regime is provided in the upper-left-hand panel of Figure 9, which is a schematic describing all regimes described in this appendix.

Fig. 9. (a)–(d) The four end-member regimes of surface meltwater response to ice-shelf tidal tilt. Left-to-right panels show the difference between freely moving water (low hydraulic resistance) and water that is highly constrained by friction. Top-to-bottom panels show the difference between a completely rigid ice shelf and an ice shelf that behaves perfectly flexibly. (e) The situation where ice-shelf rigidity is moderate (in between the two cases in the panels above) and where there is moderate resistance to surface meltwater flow. Note that dimensions on all the above illustrations are not to scale and the diagrams are schematic.

### Regime II. Friction-free, flexible ice-shelf

We now consider the same conditions as above, but with the exception ${\cal F} \ll 1$. In this case, the ice shelf is very thin (e.g. *H* < 10m), fractured (thus reducing effective plate thickness) or anelastic (e.g. Young's modulus *E* at the low end of its range, ~ ≪ 1 GPa). This regime is for an ice shelf that is much more flexible than the McMIS, and thus unlikely to exist, except in the case of sea ice (however, sea ice would unlikely be able to hold surface meltwater regimes that are kilometers in horizontal scale without draining). The equations in dimensional form are:

In this circumstance, the solution is found to be

Considering again the local linearization of the tidal regime (and choosing the *x* direction to align with the tidal tilt), we have

The effect of buoyancy is seen to amplify the meltwater-layer depth variations by at least an order of magnitude. What might thus be only an order of a centimeter or less could amplify to >37 cm. This is an important amplification that exists due to the fact that buoyancy between the ice/meltwater system and the displaced seawater below the ice shelf is in play. A schematic diagram of this behavior is shown in the lower left panel of Figure 9.

### Regime III. Diffusive meltwater boundary layer, flexible ice-shelf

In the previous two representative regimes (I and II), the resistance to water flow, as determined by the inverse of the size of the frictional damping time constant *τ*, has been assumed to be sufficiently small that the speed of adjustment to any kind of gradient on the free surface of the meltwater layer is relatively large (e.g. even with friction restraining meltwater movement, meltwater moves sufficiently fast as to be able to cross the meltwater basin several times during the diurnal or semidiurnal tidal cycle). This then gives the result that the surface-meltwater layer simply adjusts its depth so as to maintain a flat free surface even when the tide is causing the ice shelf below the meltwater layer to tilt back and forth. In the next two representative regimes (III and IV), we shall consider the consequences of large friction, small *τ*, which means that surface-meltwater will develop free-surface tilts in concert with the tide.

For representative behavior III, we assume two scaling relationships involving the way in which friction, as determined by the damping time scale *τ*, interacts with the diffusion term and forcing term when scaled according to Eqn (B1),

and additionally assume, as in Regime II, that the ice shelf is extremely flexible (${\cal F} \ll 1$). With these scaling relationships, we have made two distinctions of length scale: *L* _{m} represents the typical horizontal span of the surface-meltwater regime, and is appropriate in judging the importance of the diffusion term, and *L* _{tide} ≫ *L* _{m} represents the much larger horizontal span over which the tide elevation *T* varies, and is appropriate in judging the importance of the forcing term at points away from boundaries of the meltwater regime. The second scale relationship shown above allows us to argue that the non-homogeneous forcing term (due to the tide) can be ignored:

This, in practice, simply acknowledges that the variation in *T*(*x*, *y*, *t*) is so large in horizontal scale that only its linear terms will likely be visible in the field area. This means that the behavior in regime III is governed by a homogeneous diffusion equation:

subject to,

at boundaries, where ${\bf n}$ is the outward pointing normal to the boundary. Far from boundaries, the meltwater layer will flow in the direction determined by $\nabla T$, but this flow will have negligible divergence through most of the span of the meltwater regime, leading to ∂*η*/∂*t* ≈ 0 and *η* ≈ 0 through the interior of the regime, because *T* will be effectively a linear function of *x* and *y* having no divergence. At boundaries, however, the flow induced in the meltwater layer by the nearly uniform tilt of *T* must satisfy a no-flow boundary condition. Hence there is convergence and ∂*η*/∂*t* ≠ 0 in the vicinity of the boundary.

We can investigate the behavior of the solution near a boundary by considering one of the classic solutions to diffusion problems, one which involves a sinusoidally varying source term at the boundary of a semi-infinite domain:

subject to the boundary conditions:

where *ξ* is a local horizontal coordinate that is perpendicular to the boundary in question, taken to be positive in the direction of the interior of the meltwater regime. The boundary condition at *ξ* = 0 expressed in Eqn (B19) ensures that the free-surface elevation gradient normal to the boundary must be zero. The second boundary condition at *ξ* → ∞ represents the condition that *η* approaches the solution *η* = 0 at points far from the boundary where the lack of forcing, i.e. $\nabla ^2 T \approx 0$, ensures a homogeneous response. Physically, what is represented here is the ‘piling up’ at the boundary of a meltwater flow that is otherwise free of divergence. As the meltwater piles up at the boundary, the region where *η* ≠ 0 gradually works its way into the interior of the meltwater system in a diffusive manner.

The solution to the above boundary value problem, when $\lpar {\partial T}/{\partial \xi }\rpar \ {\rm at}\ \xi = 0$ is taken to be *T*′cos(σ*t*), is

where *σ* is the tidal frequency (e.g. ${2 \pi }/{24\, {\rm h}}$ for a diurnal tide), the *e*-folding decay scale as *ξ* increases away from the boundary is,

and the amplitude *A* is:

This idealized solution has two important properties that help to understand the meltwater regime. First, the solution encompasses a variety of phases relative to the phase of the tide, with phase varying over short distances in the direction perpendicular to the boundary. This allows us to explain time lags between the high or low tide and the timing of greatest meltwater-layer depth (maximum of *η*). Second, the solution is exponentially damped with distance from the boundary, and this provides a means of estimating an intrinsic length scale, i.e. *γ*, for surface-meltwater lake features that is independent of any pre-existing surface topography of the ice shelf. These points will be discussed further below. The two behaviors associated with Regime III are shown (for either rigid or completely flexible ice shelf) in the two right-hand-side panels of Figure 9.

### Regime IV. Diffusive boundary-layer, variably rigid elastic ice shelf

When the viscoelastic rigidity of the ice shelf is considered, as is the case when ice-shelf thickness is larger than a few meters, the role of buoyancy accommodation of the meltwater layer, i.e. the degree to which the ice shelf sinks under the weight of the meltwater layer, becomes less important, but does not vanish. The net effect is complicated, and so a simple illustrative analytic solution as derived for Regimes I to III cannot be written. We can, however, estimate an approximate solution when ice-shelf rigidity is in play by modifying the solution for Regime III to reduce the effect of buoyancy. This approximate solution differs primarily through a changed diffusivity:

where *β* ∈ [0, 1] is a coefficient that varies from 0 to 1 depending on whether the ice shelf is perfectly rigid (*β* = 0) or flexible (*β* = 1), respectively. The only significant change to the idealized behavior of the meltwater system from the situation in Regime III (which implicitly has *β* = 1) is that the exponential-decay scale is increased, and the amplitude is decreased, as *β* decreases toward 0 (perfect ice-shelf rigidity):

and,

Another complexity of the behavior associated with the variably rigid ice shelf is that the domain of non-zero deflexion of the ice-shelf neutral plane, ζ ≠ 0 extends beyond the area covered by the surface meltwater layer. This allows the possibility that two or more separate, disconnected meltwater regimes could influence each other during a tidal cycle.

To illustrate the effects of variable ice-shelf rigidity over a range of parameter values, we resort to using a numerical solution (implemented using the finite-element package COMSOL) to couple the ice-shelf viscoelastic flexure to the meltwater layer in a 1-dimensional geometry consisting of a horizontal span of meltwater sufficiently large as to isolate the diffusive boundary layers discussed in the previous section. Figure 10 shows several representative solutions to the coupled equations simulated using

or about 1 cm of tidal tilt per kilometer (consistent with local conditions on the Larsen B Ice Shelf), *σ* _{M2} = 2π/12.42 h^{−1} is the tidal frequency for the dominant semidiurnal tide in the Weddell Sea (Robertson and others, 1985; King and others, 2011), *τ* = 500 s, which is a fairly strong frictional damping term associated with a Manning roughness parameter of about 0.15 m s^{1/3}, and the ice-shelf is treated with a purely elastic ice-shelf rheology (viscous behavior will be examined later) with a flexural rigidity corresponding to various ice thicknesses (10, 50, 250 m) and a Young's modulus *E* of 5 GPa.

Fig. 10. Meltwater-layer depth perturbation (*η*, cm, blue line), ice-shelf deflexion due to purely elastic flexure (*ζ*, cm, red line), tide (*T*, cm, dashed black line), free-surface of meltwater layer (blue dots) and ice-shelf stress at the upper surface (*T* _{xx}, kPa, black line) for three representative ice thicknesses (*H* = 10,50 and 250 m, panels (a–c), respectively). The time for which the solution is shown is π/4 units of phase after low tide (when the slope of *T* is maximum), and is chosen to show the depth of the meltwater-layer perturbation at its largest. Note the change in scale for each panel.

One of the key features of the elastic flexure of the ice shelf in response to the tidally forced meltwater fluctuations is the creation of flexure stresses which, as a result of the thin-plate approximation, are strongest, and of opposite sign, at the upper and lower ice-shelf surfaces (the upper and lower skin surface) and which vary linearly with vertical distance in the ice shelf. Figure 11 displays the maximum tensile stress attained at the upper surface of the ice shelf as a function of ice thickness (varying from 5 to 500 m) and for a variety of rheological treatments. The five rheological treatments are: (1–3) pure elastic, with *E* = 1,5,10 GPa, representing the range of natural ice, and (4–5) viscoelastic (following the method of MacAyeal and others, 2015), with *E* = 5 GPa and ν = 1 × 10^{13} and 2.2 × 10^{14} corresponding to Maxwell times of about 30 min and 12.42 h (the M2 period), respectively.) The maximum tensile stress *T* _{xx} achieved in the 1-D numerical model used to derive the results in Figure 10 increases with increasing Young's modulus (*E*) but decreases with *H*, as is displayed by the rightward decay of the curves in Figure 11. The flexural rigidity in the purely elastic treatment, *D*, is proportional to *EH* ^{3}; so the decay with increasing ice thickness reflects the fact that the length scale over which non-zero bending moments are induced in the ice shelf by the shifting meltwater load is an increasing function of *H*. Only oneviscoelastic example is shown, and is produced for an unrealistically low ice viscosity (10^{13} Pa s) associated with an unrealistically short Maxwell time ν/*E* of approximately 30 min. For Maxwell times at or above the diurnal and semidiurnal tidal periods, the viscoelastic behavior is indistinguishable from the purely elastic behavior.

Fig. 11. Maximum tensile stress achieved at the time of greatest meltwater-layer depth as a function of ice thickness *H* and for various values of the Young's modulus *E* under the assumption of pure elastic rheology. Two viscoelastic solutions are also plotted. One with *E* = 5 GPa and a Maxwell time of 12.42 h (the periodicity of M2) plots exactly above the line for the purely elastic case with *E* = 5 GPa. The second with *E* = 5 GPa and assuming an ice viscosity of 1 × 10^{13} Pa s (an extremely low value, corresponding to a Maxwell time of 0.55 h) is located between the two elastic cases of *E* = 1,5 GPa. The maximum tensile stress is achieved at the surface of the ice shelf at the boundary of the meltwater regime near *x* = 0.

The magnitude of the maximum tensile flexure stress of the ice shelf may exceed what may be capable of producing fracture in the skin (upper surface) of the ice shelf (Banwell and others, 2013). This is particularly significant given the fact that this tensile stress is cyclic and thus giving rise to the possibility of fatigue failure which occurs at stress levels that are reduced below the fracture criterion for time-independent stresses. Although the 1-D results shown in Figure 11 are idealized both in terms of geometry and tidal forcing (e.g. the meltwater-layer extends infinitely in the direction perpendicular to *x*, and the tidal-tilt magnitude d*T*/d*x* = 10^{−5} is at the high-end of what is feasible in natural settings), ice shelves with an effective elastic plate thickness of ~50 m or less will experience ~50 kPa tensile stresses (in their skin) repeatedly with each tidal cycle. This suggests that tidally induced flexure stresses due to movement of surface meltwater will have a bearing on the stability of thin ice shelves such as the McMIS (which ranges between 20 and 100 m in thickness, Rack and others, 2013). Thick ice shelves, even ones with extensive surface meltwater systems, such as the Larsen B Ice Shelf prior to its 2002 collapse, will experience maximum tensile stresses on the order of 5–20 kPa, which are roughly 5–20% of the single-cycle fracture criterion of ~100 kPa. Cyclic stresses in the 10% range, after sufficient number of cycles, could be sufficient to fracture an ice shelf as is consistent with the basic principles of fatigue failure.

## APPENDIX C. Estimation of tidal tilt on the McMurdo ice shelf

Over ice-shelf covered seas at locations distant from coasts and grounding lines, dry firn-covered ice shelves ride passively atop the ‘ocean-scale’ barotropic tidal waves that make up the semidiurnal, diurnal and long-period response of the global ocean to the astronomical tide potential (Platzman, 1978; Henderschott, 1980). Counterintuitively, the great thickness of the ice shelf, sometimes far thicker than the ocean layer which is trapped below it, does little to influence the propagation of the ocean tide beyond determining the effective depth of the ocean (actual seabed depth minus ice-shelf draft) and providing a second rigid surface at which to generate friction (along with the sea bed). This characteristic results from the fact that long-period barotropic ocean tides propagate as shallow water waves (MacAyeal, 1984; Padman and others, 2003b, 2018) for which vertical inertia of the ice and ocean beneath exerts a negligible influence on pressure (Pedlosky, 1987; Hutter and others, 2014). For shallow water waves, only horizontal components of inertia are relevant, and horizontal motions of ice shelves driven by tide are easily disregarded when considering tidal wave propagation (Brunt and DR, 2014).

Near coasts and grounding lines, however, the effect of the ocean-scale tide on vertical ice-shelf motion is complicated by at least four factors: (1) elastic rigidity of the ice shelf leading to flexural boundary layers within ~10 km of grounding lines (e.g. Vaughan, 1995; Rosier and others, 2017), (2) decreasing water depth and increasing friction as the water column trapped between the sea bed and the ice-shelf narrows (Georgas, 2012; Rosier and others, 2014), (3) the effects of nonlinear inertia terms that create asymmetry and generate harmonics (Aubrey and Speer, 1985; King and others, 2011), and (4) local resonance (Arbic and others, 2009). These effects generate perturbations from the large-scale character of the ocean-scale tidal wave propagation which are typically restricted to relatively small coastal regions. The rapidly increasing thickness of the ice upstream of a grounding line inhibits the rise and fall of the ice shelf near the grounding line, which leads to a kilometer-scale boundary layer along the grounding line where the ocean-scale tide is damped to zero by viscoelastic forces in the ice (e.g. Fricker and Padman, 2006; Walker and others, 2013; Rosier and others, 2017).

Ice-shelf tilt is defined by $\nabla \zeta$, where *ζ* is the perturbation of the neutral surface of the ice shelf (typically half way between the surface and bottom) about which elastic and viscoelastic deformation occurs, and where ∇ is the gradient with respect to horizontal coordinates. This tilt has a small component caused by the ocean-scale tidal wave in regions away from coasts and a large component caused by processes operating in grounding zones and where water-column thickness is restricted. Ocean-basin scale tilt is caused by two factors: (1) the variation of tidal phase in the direction of tidal-wave propagation and is typically parallel to the basin-scale trend of the ocean coastline (e.g. Platzman, 1978; Henderschott, 1980), and (2) variation of tidal amplitude in a direction perpendicular to the coast caused by both coastal amplification due to shallowing and the effects of Earth's rotation (Pedlosky, 1987). The Kelvin wave (Platzman and others, 1981; Pedlosky, 1987) is a good idealized model of what the ocean tide looks like. These two variations are seen, for example, in the dominant semidiurnal M2 tide's amplitude and phase maps for the Southern Weddell Sea (Robertson and others, 1985; Padman and others, 2003a, 2018; King and others, 2011) and in the dominant diurnal K1 tide's amplitude and phase maps for the entire Southern Ocean (Platzman and others, 1981; Padman and others, 2018). The wavelength of the large-scale, ocean-scale tide is typically >6000 km; thus, for a tidal amplitude of 1 m, the maximum tilt in the direction of propagation is about 10^{−7}, or about 1 cm in 100 km. Maximum tilt in the direction perpendicular to phase propagation is determined by the Rossby radius of deformation (Pedlosky, 1987), ${\sqrt {gh}}/{f}$, where *g* is the acceleration of gravity at the Earth's surface, *h* is mean water-column thickness (depth in open water) and *f* is the magnitude of the Coriolis parameter for the latitude in question, as well as by small-scale energy continuity considerations associated with varying *h*. For the Weddell Sea, the typical scale of tide variation perpendicular to the coastline is about 1 cm in 1–10 km (Robertson and others, 1985; Padman and others, 2003a; Mueller and others, 2012; Zaron, 2018). This implies a tilt of about 10^{−5} to 10^{−6}.

In the kilometer-scale boundary layer associated with grounding zones, ice-shelf tilt is on the order of 10's of cm per kilometer or 10^{−4} (e.g. Rosier and others, 2017). The grounding zone thus corresponds to a location where surface slopes can be as large as for ablation zones at the edges of the Greenland ice sheet. The tilt in the grounding zone is also the largest sea-surface tilt in the World Ocean, exceeding, for example the slope of dynamic topography, ~ 10^{−5}, at the edges of the Gulf Stream off the Eastern coast of North America (Knudsen and others, 2011; Rosier and others, 2017). While little has been done to observe tide asymmetry and harmonic excitation in ice-shelf-covered seas, observations of such phenomena elsewhere suggest that these effects could lead to centimeter-scale variation of ice-shelf surface elevation over kilometer-scale horizontal distances, thus also about ~ 10^{−5} (e.g. Aubrey and Speer, 1985).

On discovering the diurnal variability in meltwater depth at Rift Tip Lake, the question we considered first was whether the cause could be the large-scale ocean tide in the Ross Sea, which is primarily diurnal (MacAyeal, 1984; Padman and others, 2003a, 2003b). To investigate whether an ice-shelf surface tilt was associated with the tide propagating through McMurdo Sound, we examined GPS data we collected during the 2016/17 field season when the GPS array covered a larger area. GPS data collected in 2015/16 were from an array that was horizontally confined (within 1 km) and thus unable to allow a tilt signal to emerge from the GPS noise (of several centimeters). (We use the 2016/17 data because the location of GPS receivers in the second field season was more optimal for detecting tidal amplitude variation, and because the tidal propagation in McMurdo Sound is likely to not vary significantly from year to year.) Two GPS stations (labeled 1 and 2 in Fig. 2B, and identified as Ring GPS3 and Peanut GPS3 in Banwell and others, 2019) are also located on the relict superimposed lake ice that now forms pedestals that rise above the ice-shelf surface due to differential ablation (Macdonald and others, 2019), which is where we expect the GPS data to be least influenced by water movement during the early melt season (i.e. we expect these stations to reflect the tidal motion of the ice shelf without the influence of flexure in response to water loads). Description of the GPS measurements, data processing and uncertainty is provided in Banwell and others (2019).

Figure 12a shows the GPS height above geoid (mean value subtracted and smoothed with a 3-h running mean) and (bottom panel) the difference between GPS1 and GPS2 as a function of the height above geoid at GPS1. The result of fitting a line through the noisy GPS data in Figure 12b suggests that the tide amplitude differs by about 2 cm between Ring Lake (greater amplitude) and Peanut Lake (lesser amplitude). These two stations are about 3 km apart. This variation in tidal amplitude forms the basis for our use of 0.5 (≈2/3) cm per km as the tidal tilt in the vicinity of Rift Tip Lake for the exploratory numerical experiments conducted in this study.

Fig. 12. (a) Comparison of GPS-observed elevation over a 4-d period (13–17 January 2017) between Ring and Peanut lake GPSs (GPS 1 and 2 in Fig. 2b, also named Ring GPS3 and Peanut GPS3 in Fig. 1 of Banwell and others, 2019). (b) Difference between GPS elevations at Ring and Peanut against tidal elevation at Ring over the 4-d period. Blue line indicates zero difference, red line indicates a least squares fit through all data points suggesting that the tidal amplitude is approximately 2 cm greater at Ring than at Peanut. The two GPS stations are 3.1 km apart. This forms an observational basis for our use of 0.5 (≈2/3) cm/km as the tidal tilt in the numerical experiments of this study.