The discontinuous Galerkin model of shallow water tides

The DG model of shallow water tides was first described by Salehipour et al. (Reference Salehipour, Stuhne and Peltier2013) where it was benchmarked against the aquaplanet test cases 2 and 5 described by Williamson et al. (Reference Williamson, Drake, Hack, Jakob and Swarztrauber1992) for shallow water models as well as via the simulation of the evolution of a baratropically unstable jet due to perturbation of an initially stable zonal flow, as was originally proposed as an additional shallow water theory benchmark by Galewsky et al. (Reference Galewsky, Scott and Polvani2004). The DG model has also been validated against the TOPEX/Poseidon data for the M2 and K1 tidal constituents under modern conditions, and the results for the M2 and K1 tidal constituents at last glacial maximum (LGM) have been presented by Griffiths and Peltier (Reference Griffiths and Peltier2008, Reference Griffiths and Peltier2009). The methodology, briefly summarized, involves using the 2D spherical shallow-water tidal equations (Egbert et al., Reference Egbert, Ray and Bills2004) in their flux form as:

(1)$$\partial _{t}{\bf Q} + \nabla \cdot {\bf F}( {\bf Q} ) = {\bf S}( {\bf Q} ) .$$

Here Q represents all conserved quantities, F = [F^x, F^y, F^z] is the flux term, S is the source term, and Q = [φ φ u]^T,

(2)$${\rm F}^{\rm x} = \left[\matrix{\phi u \hfill \cr \phi u^2 + \displaystyle{1 \over 2}\phi^2 \hfill \cr \phi uv \hfill \cr \phi uw \hfill} \right], \;{\rm F}^{\rm y} = \left[\matrix{\phi v \hfill \cr \phi uy \hfill \cr \phi v^2 + \displaystyle{1 \over 2}\phi^2 \hfill \cr \phi vw \hfill} \right], \;{\rm F}^{\rm z} = \left[\matrix{\phi w \hfill \cr \phi uw \hfill \cr \phi vw \hfill \cr \phi w^2 + \displaystyle{1 \over 2}\phi^2 \hfill} \right]$$

and

(3)$${\rm S} = \left[\matrix{{0, \;-f\,{\hat{{\rm e}}}_g \times ( \phi {\rm u}) + \phi \nabla ( \phi_{eq} + \phi_{sal}-\phi_s\rpar} {-\displaystyle{g \over \rho}( D_{BL}+D_{IT}\rpar}} \right]^T$$

Here φ = gh(x, t) is the geopotential height, ρ is the density, g is the gravitational acceleration, φu is the depth-averaged horizontal transport, u = (u, v, w)^T is the velocity vector, ê _g is a unit vector pointing normal to the surface of the sphere, f is the Coriolis parameter, φ_eq is the equilibrium tidal variation induced by astronomical bodies, φ_sal is the influence on the tide due to gravitational self-attraction and loading, and φ_s is the equilibrium geopotential height. The fluid depth, h(x, t), can be expressed as the sum of the depth measured from the undisturbed ocean surface to the sea floor h_s(x, t) and free-surface elevation ζ(x, t). As the velocity vector is expressed in Cartesian coordinates, it must be explicitly constrained to the sphere using the same method based on the use of projection matrices as was employed by Williamson et al. (Reference Williamson, Drake, Hack, Jakob and Swarztrauber1992). D_IT and D_BL are the dissipation due to internal tides and the ocean floor boundary layer, respectively.

The form of the tidal variation due to the gravitational attraction of distinct astronomical bodies, the influence on the tide due to self-attraction and loading (Accad and Pekeris, Reference Accad and Pekeris1978), as well as boundary layer (Taylor, Reference Taylor1919) and internal tidal dissipation (Griffiths and Grimshaw, Reference Griffiths and Grimshaw2007) used in this investigation are identical to those used by Velay-Vitow et al. (Reference Velay-Vitow, Peltier and Stuhne2020) following the work by Salehipour et al. (Reference Salehipour, Stuhne and Peltier2013). They are given respectively by:

(4)$$\eqalign{\phi _{eq} = \,& g\sum\limits_l {A_l {\kappa}_l\left({\displaystyle{1 \over 2}-\displaystyle{3 \over 2}{\sin }^2( \theta ) } \right)} \cos ( \omega _lt) \cr & + g\sum\limits_d {A_d {\kappa}_d\sin ( 2\theta ) } \times \cos ( \omega _dt + \lambda ) \cr & + g\sum\limits_s {A_s {\kappa}_s{\cos }^2( \theta ) } \cos ( \omega _st + 2\lambda ) , \;} $$

(5)$$\phi _{sal} = ( 0.085g) \zeta $$

(6)$$D_{BL} = \rho c_d\vert u\vert u, \;\,c_d = 0.0025$$

and

(7)$$D_{IT} = \displaystyle{{\rho {\bar{N}}^2h_s} \over {3\omega }}( u\cdot \nabla h_s) \nabla h_s \times \left\{{\matrix{ {1, \;} & {\vert f\vert \,< \,\omega } \cr {0, \;} & {\vert f\vert \,> \,\omega } \cr } } \right.\cdot $$

Here θ and λ are latitude and longitude, A is the equilibrium amplitude, ω is the tidal frequency, and κ is the appropriate surface load Love number (Farrell, Reference Farrell1972). These values can be found in table 1 of Arbic et al. (Reference Arbic, Garner, Hallberg and Simmons2004). Finally, N̄̄² is the vertically averaged buoyancy frequency of the ocean. For the purposes of this analysis, the longitude- and latitude- dependent vertical average of the buoyancy frequency of the oceans is assumed to be identical to that for LGM obtained on the basis of a coupled climate model simulation of glacial conditions (Vettoretti and Peltier, Reference Vettoretti and Peltier2013). Though the question of deep stratification has not been completely resolved in the literature, the DG model employed is only very weakly dependent on the buoyancy frequency. As such, the exact choice of how to represent these frequencies does not affect our results.

Table 1.

Kinetic and potential energies of the M2 tidal constituent prior to and after Heinrich event 3 (H3).

Review of Velay-Vitow et al. (Reference Velay-Vitow, Peltier and Stuhne2020)

Velay-Vitow et al. (Reference Velay-Vitow, Peltier and Stuhne2020) advance a model of tidally triggered marine ice stream instability as a possible explanation for H-events. In this model, the strait is initially assumed to be filled with an ice stream that is responsible for draining a significant region of the eastern Laurentide Ice Sheet. Additionally, the M2 tidal constituent is in a strongly enhanced state of resonance in the North Atlantic. Constructive interference of the M2 and S2 tidal constituents (the spring-neap cycle) is suggested to be responsible for triggering an ice stream instability in the Hudson Strait, which proceeds rapidly inland until it is arrested at Southampton Island. The underlying crust, depressed under the weight of the ice stream, must be of sufficiently steep retrograde slope for this instability to propagate (Schoof, Reference Schoof2012). With the weight of the ice removed, the floor of Hudson Strait and the surrounding region begins to rebound under the action of the isostatic adjustment process (e.g., Peltier et al., Reference Peltier, Argus and Drummond2015, for recent review) until such a time as the slope is sufficiently reduced to allow a re-advance of the Hudson Strait ice stream to occur. This ice stream, no longer prone to tidal destabilization, remains in place and continues to thicken until the weight of the reestablished ice stream is able to depress the Earth's crust sufficiently to render the ice stream susceptible to a further deglaciation event. It is in this fashion that the authors imagine the approximately 7000-year pseudo-periodicity of H-events to be achieved. This process is illustrated by the cartoon in Figure 2. The proposed model is suggested to be an example of Self Organized Criticality (SOC), a phenomenon in which non-linear systems may be characterized by the existence of a point of instability in the space of its parameters, a point which is also an “attractor” of its dynamics. That is to say, non-linear systems are continuously drawn toward a state in which they experience a rapid change in behavior.

Figure 2.

(color online) In panel one, the strait is shown prior to a Heinrich event; the underlying crust slopes up to the Atlantic, and the strait is filled with grounded ice. The North Atlantic resonance is strong, and the amplitude of the tides is large. In panel two, a Heinrich event has occurred, and icebergs laden with ice-rafted debris are released into the Atlantic Ocean. As the ocean is no longer in a state of maximum tidal resonance, the amplitude of the tides is decreased. Panel three shows the crust rebounding as it is no longer depressed under the weight of the ice. In panel four, the ice begins to reform as the slope of the underlying crust is no longer sufficiently steep to prevent ice stream readvance. Note that the floating ice stream that could potentially exist at the margin of the ice stream is not shown. Figure and caption adapted from Velay-Vitow et al. (Reference Velay-Vitow, Peltier and Stuhne2020).

This concept was introduced into the physics literature by Bak and Tang (Reference Bak and Tang1989) who discussed it in terms of a “sand pile” analogy. In this thought experiment sand grains are continuously dropped onto a horizontal surface. Eventually a conical “pile” of sand is created in which the surface slope is at the “critical angle of repose.” As further grains of sand continue to fall onto the pile, “avalanches” of various size eventually occur at points on the conical surface that are forced to exceed the critical angle of repose. These act so as to restore the local slope of the pile back to the critical angle. This process simply repeats. Avalanches are certainly not periodic but rather episodic, and of variable strength in this particular model. During MIS 3 the SOC terminology appears to be well motivated. Individual Heinrich events are analogous to the avalanches of the sand pile (the eastern sector of the Laurentide Ice Sheet). The continuous forcing of the sand pile by the falling grains of sand is analogous to the tendency of the Hudson Strait ice stream to reform following the occurrence of each Heinrich event instability; the straight is continually “forced” by the flow of ice into it. This continuous reformation of the ice stream is supported by ice dynamical processes that cause ice to flow back into the void created by instability of the ice stream once the slope of the bed has been sufficiently reduced by the glacial isostatic adjustment process.

Salehipour et al. (Reference Salehipour, Peltier and Caulfield2018) have provided a recent fluid mechanical example of SOC. These authors further reference what appears to have been the first such example in the geophysical literature, namely that discussed in Butler and Peltier (Reference Butler and Peltier1997) in their analysis of the problem of convective mixing in the deep earth in which the presence of an endothermic phase transition enforces an intense time dependence upon the mixing process. Although in many examples of SOC there is a lack of a fixed time scale for recurrence, the controlling factor we are positing is the depression and uplift of the Earth's crust, which is expected to impress an approximate characteristic repeat time on the recurrence of the instability.

An important additional result of the analysis presented by Velay-Vitow et al. (Reference Velay-Vitow, Peltier and Stuhne2020) was that the removal of the ice grounded in Hudson Strait had a significant effect on the state of the North Atlantic resonance of the M2 tidal constituent as shown in Figure 3, in which we show the difference between the amplitude of the M2 tide before and after instability occurrence. If it is imagined that an altogether different ice stream located in the ice-covered continental region of the eastern sector of the North Atlantic was responsible for anomalous HEs, as was suggested in Hemming (Reference Hemming2004), it would be helpful to identify candidate ice streams and to consider whether or not they might be susceptible to the same instability mechanism as the one we have shown to be highly plausible for the majority of H-events. If it were found to be the case that eastern-sector instabilities were highly plausible, this would undermine the plausibility of our suggested mechanism for H-event recurrence time. Our goal in this paper is to provide detailed examinations of several candidate ice streams and to comment on the question as to whether they should be viewed as plausible possibilities for the explanation of anomalous H-events.

Figure 3.

In our model, the effect of H1 on the resonant state of the M2 tide in the North Atlantic was particularly profound outboard of the Hudson Strait, with a drop in tidal amplitude of between 0.3 and 0.6 m. Figure and caption adapted from Velay-Vitow et al. (Reference Velay-Vitow, Peltier and Stuhne2020).

Candidate models for Heinrich event 3

In modeling H3, the initiation and duration of the instability and the extent of the deglaciation must be determined. We make use of the results in Hemming (Reference Hemming2004) in the design of a plausible model and choose 31 ka as an initiation time for the (assumed to be related) ice stream instability and 1200 years as the duration of the model event. The duration of the modeled event was determined by averaging the three different estimates for the duration of H3 listed in table 3 of Hemming (Reference Hemming2004). To determine the plausible regions of ice stream instability that might have been the source of a non-Hudson Strait H3 event, we first examined the crustal depression predicted at each of the most plausible candidate locations by the ICE-7G NA (VM7) model of Roy and Peltier (Reference Roy and Peltier2017, Reference Roy and Peltier2018) to search for similar features to those of the Hudson Strait reported by Velay-Vitow et al. (Reference Velay-Vitow, Peltier and Stuhne2020). Specifically, we search for regions in which the nature of the crustal depression beneath grounded ice cover is such that the depression increases with distance “up ice” from the point at which grounded ice cover abuts the Atlantic Ocean. The nature of the crustal depression beneath grounded ice in the ICE-7G NA (VM7) model of Roy and Peltier (Reference Roy and Peltier2017, Reference Roy and Peltier2018) for each of our candidate locations is presented in Figure 4. Panel A of the figure illustrates the depth of the crustal depression in Hudson Strait, and the remaining panels of the figure show the same data for our three eastern North Atlantic candidate ice streams. This characteristic of the crustal depression beneath an ice stream has been demonstrated by Schoof (Reference Schoof2007, Reference Schoof2012) to be the most important requirement for an ice stream to be susceptible to instability when the grounding line is subjected to a sufficiently intense buoyancy forcing to move it downslope. In all three of the eastern North Atlantic candidate ice streams, the crustal depression along the bed of the proposed ice stream is smaller in magnitude than that observed in the Hudson Strait case, suggesting that the requirement for instability will be more difficult to meet in these regions.

Figure 4.

The crustal depression under the weight of terrestrial ice sheets at each of the three candidate eastern North Atlantic locations, as well as the Hudson Strait, immediately before the onset of H3 at 31 ka is presented. In panel A), the crust underlying the Hudson Strait has qualitative features that are nearly identical to those identified by Velay-Vitow et al. (Reference Velay-Vitow, Peltier and Stuhne2020). The Irish Sea crustal depression in panel B) shows only a slight inland deepening, while the Norwegian ice stream, presented in panel C), shows a clear increase in depth in the up-strait direction. As shown in panel D), the center of the Barents Sea is depressed relative to the edges, and the depth increases from west to east. Furthermore, there is a local maxima of depth separating the North Atlantic Ocean from the Barents Sea. All three candidate locations have features somewhat similar to those observed in the Hudson Strait. The black bars indicate the approximate location of candidate ice streams.

Velay-Vitow et al. (Reference Velay-Vitow, Peltier and Stuhne2020) suggested that the forcing of the ice stream at its point of contact with the ocean, which was needed to induce instability, might have been provided by anomalously high tides, a condition that was extremely well satisfied by the conditions of the Hudson Strait ice stream prior to H1. These conditions are qualitatively similar to the conditions in the Hudson Strait prior to H3 (see Figure 4). If similar features were found to exist at one or more of our eastern Atlantic candidate locations, this would allow the same model employed in the H1 analysis to be applied to these additional candidate ice streams. In addition to the three locations on the east coast of the North Atlantic, the Hudson Strait ice stream was also examined. In each case, H3 was modeled by modifying the ICE-7G NA (VM7) model of ice-loading history from Roy and Peltier (Reference Roy and Peltier2017, Reference Roy and Peltier2018) to determine the paleobathymetric depth that would exist once the grounded ice was eliminated from the region assumed to have been filled by the candidate ice stream; this depth is then used as an input to the shallow-water tidal solver. Although this model includes no ice dynamical processes, such as calving, the grids generated for the tidal solver only depend on whether a given cell is ocean or land. The additional sea-level rise associated with realistic ice stream dynamics is assumed to be negligible and is therefore ignored.

As in Velay-Vitow et al. (Reference Velay-Vitow, Peltier and Stuhne2020), we make use of a non-uniform grid that is more highly resolved in each of the proposed locations, as well as in the region of high resonance in the Labrador Sea, to accurately model the effect deglaciating the proposed ice streams has on the tidal regime in the North Atlantic. The grids we employ have a maximum internodal spacing of 30 km in the open ocean, and a minimum of 7.5 km near the coasts. The grid used for the Ireland case is presented in Figure 5.

Figure 5.

(color online) The globally unstructured grid is locally refined in the vicinity of the hypothesized ice stream in the Irish Sea. In the upper left of the figure, the transition region from comparatively higher internodal spacing of 30 km to the lower spacing of 7.5 km is seen.