Co-seismic eruption and intermittent turbulence of a subglacial discharge plume revealed by continuous subsurface observations in Greenland

In the Arctic, subglacial discharge plumes have been recently recognised as a key driver of fjord-scale circulation. However, owing to the danger that accompanies prolonged observations at plumes, no time-series data are available. Here, we present results showing the chaotic and irregular dynamics of a plume revealed by continuous subsurface monitoring directly on the calving front of a Greenlandic glacier. We found intense fluctuations in the current and scalars (temperature and salinity), recognised shallow and deep tidal modulation and anomalies due to co-seismic drainage of an ice-dammed lake via the plume, and observed rapid and marked changes in stratification. Our analysis uncovers energy cascade intermittency with coherent structures, corresponding to upwelling pulses of warm water. Prior to our research, in situ evidence of time-variable plume dynamics was absent and limited to snapshots, therefore, our study and approach will enable researchers to transition from an episodic view of a plume to a continuously updated image. Direct observations of a glacial plume on the calving front reveal intense, tide-modulated fluctuations in physical properties of water, drainage of an ice-dammed lake, sudden stratification shifts, as well as coherent upwelling pulses of warm water.

P ower-plant chimneys, sewage pipes, forest fires, volcanoes and hydrothermal vents are examples of artificial and natural sources of buoyancy-driven turbulent convective plumes (also referred to as gravitational convection 1 ). At the fronts of marine-terminating glaciers, plumes are produced by a discharge of low density, buoyant meltwater from conduits at the glacier base. The density contrast between fjord-bottom saline water and fresh subglacial discharge continuously drives buoyancy at the submerged foot of a glacier 2 . The structure of the plumes depends on discharge rates, ice-cliff/conduit geometry and stratification [2][3][4] . A growing body of literature indicates a significant influence of upwelling discharge plumes on the dynamics of proglacial fjords 2,5-9 , which are becoming the last refuge for some Arctic animals as sea-ice continues to decline and glaciers retreat to land [10][11][12] . Rising plumes entrain nutrient-rich seawater 9 , and then melt and undercut the glacier front causing calving 7,8 , the mixing and upward transport of heat, nutrients and zooplankton 12 , and driving vertical and horizontal fjord-scale circulation 2,5 . Some studies suggest that when fresh water meets deep saline water, zooplankton may experience a so-called 'osmotic shock' 10,12 , while fish can be entrained by convection 13 . Their consequent rapid vertical ascent to the surface makes them easy prey, and leads to the formation of foraging hotspots for seabirds and seals 10,11,14 . Interestingly, subglacial meltwater plumes might serve as a transport mechanism for seals, who apparently use the plumes for returning to the surface 14 . The reasons attracting other endemic Arctic marine mammals, such as narwhals, to calving fronts have not yet been established 15 .
Numerical modelling of plumes and the documenting of water properties by unmanned vehicles [16][17][18] , helicopters 4,16,19 and biologging 14 have clarified some of these processes. Observational and numerical work, which aimed to estimate the response of glacier systems to submarine melt, have indicated the necessity of taking direct measurements in the plumes 2,3 . However, all efforts thus far have given only instantaneous, non-continuous data 4,5,14,[16][17][18][19][20] . To our knowledge, continuous, in situ monitoring of the plumes with a sufficiently high temporal resolution, allowing characterisation of abrupt and turbulent dynamics, is lacking. This observational gap is due to the extreme difficulty and danger in conducting continuous fjord observations near the glacier terminus. The crevassed terminus is regularly losing mass through catastrophic iceberg calving events, and the nearby sea surface may be covered by an ice mélange. Here, we resolve this disparity using comprehensive high-temporal-resolution observations for July 2017. We deployed sensors from the calving front of Bowdoin Glacier (Kangerluarsuup Sermia in Greenlandic language) directly into the plume surface footprint (Fig. 1). This 3-km-wide glacier in Northwest Greenland (Fig. 1) is a typical marine-terminating glacier, grounded in a 250-m-deep fjord with a tide-modulated ice speed of 1-3 m d −1 in summer 21,22 . The operation was conducted between two major, kilometer-scale calving events on 8 July and 1 August (Fig. 1c-e). We hung two sets of sensors from a 30-m-high ice cliff and connected them with individual cables to a recording system that was ice-screwed to the glacier surface (Fig. 1c). Every 10 s for 12 days (13-24 July) the first, 'deep' sensor recorded temperature, conductivity and pressure at a depth of~100 m, and the second, 'shallow' sensor logged temperature and conductivity at a depth of~5 m. On July 14, we directly observed, and our high-frequency time-lapse cameras, seismometer and lake-water-pressure sensor recorded, the subglacial drainage of an ice-dammed lake (hereafter Glacial Lake Outburst Flood, GLOF) via the plume with immediate impacts on its intensity (Fig. 2). The lake was located~2 km from the plume, with an area of 27,000 m 2 (Fig. 1b) and a water temperature of +1.9 ∘ C.
Understanding how plumes behave and evolve is one of the most compelling challenges in present ice-ocean interaction studies, and here we perform a synergistic analysis of this first-ofits-kind dataset. For our assessment, we acknowledge the spiky (intermittent) appearance of the time-series, which indicates turbulence and could mistakenly be interpreted as being stochastic rather than deterministic ( Supplementary Fig. 1). Therefore, we employed a combination of statistical, nonlinear, wavelet-based, spectral, seismic and other signal/image-processing approaches, as detailed in Methods, together with a description of our observations and data handling.

Results and discussion
Continuous plume monitoring: GLOF anomaly, an abrupt stratification change and tidal modulation. The outburst of the same ice-dammed marginal lake was witnessed on 6 July 2016 18 , supported by a photogrammetrically inferred lake volume (~190,000 m 3 ) and a computed mean daily melt rate of glaciers discharging meltwater into Bowdoin Fjord through three plumes (~6.4 × 10 6 m 3 d −1 or~265,000 m 3 h −1 or~74 m 3 s −1 ) 18 . In July 2017, the lake was approximately twice its former size due to differences in the shapes of the ice margin and the lake, and only two plumes were visible (Fig. 1b). Assumption of a drainage of 380,000 m 3 (i.e. a twice larger volume than in 2016) within 8.5 h at a constant rate yields a subglacial discharge of~45,000 m 3 h −1 (or~12 m 3 s −1 ). The partitioning of total meltwater discharge between the two plumes is unknown. A 1:1 partitioning would imply that the flux from the GLOF could increase the average subglacial discharge driven by surface melting 18 by an additional flux of~34%. This value is non-negligible, and could suddenly overwhelm the subglacial drainage system, especially as it might be considered the lower boundary. More specifically, (1) the volume of the lake could be larger than the assumed one (the depth of the lake increases on the ice side due to the valley slope); (2) the baseline drainage rates via the monitored plume could be smaller than 50% of the total because it had a smaller surface footprint than the second plume (Fig. 1b); and (3) the peak discharge rate from the lake could be non-constant and higher than the estimated one. For example, the seismic signal of the GLOF increased with time (Fig. 2c). Moreover, time-lapse imagery shows that concentric surface ripples from the plume intensified during the GLOF (after 22:40, 14 July; Supplementary Movie 1).
Trajectories of the deep sensor in temperature-salinity-pressure space (Fig. 3) show recurrent dynamics mixed with highly irregular motion triggered by the GLOF with several ejections and departures of the sensor towards the surface by iceberg pull and current pulses (Supplementary Movies 2-4).
Overall, the deep sensor exhibited a narrow variance in temperature and salinity, and showed a drop in salinity of 1 PSU lasting at least half a day after the GLOF (Fig. 3). For more than 80% of the time-series duration, the sensor probed depths of >80 m where we observed temperature oscillation between −1.5 and −0.5 ∘ C and an average salinity was 33.66 PSU (Fig. 3d-f). As the GLOF-induced outward push of the plume was halted, water and floating ice rushed back to the calving front (in accordance with a plume-subduction theory 23 ), causing surface warming and a freshening for 3 h (Fig. 3a- The shallow sensor showed wide variation in temperature and salinity ( Fig. 3d and e), with a notable +1.8 ∘ C change on 16 July and freshening afterwards ( Fig. 3a and b). This coincided with a 2 ∘ C increase in the mean air temperature and arrival of new fjordsurface water to the area, which replaced ice mélange with outside sea-ice floes ( Fig. 4 and Supplementary Movie 1).
Surprisingly, this surface transition was associated with a dramatic change in the motion dynamics of the deep sensor. In particular, we detected a shift in the statistical properties of the rate of pressure change (e.g. positive for descent, negative for ascent; Supplementary Fig. 1a). To reveal this feature, we first removed all amplitude information completely (i.e. data spikes) by applying one-bit normalisation in which signal amplitudes became ±1, which is a common technique in seismology. Second, we focused on the temporal variability of such statistical measures as skewness, representing the asymmetry of variable distributions.
One-bit-normalised skewness for the rate of pressure change ( Supplementary Fig. 1c) was notably negative after 16 July, implying that the sensor spent more time ascending amidst rare periods of rapid descent. Before this date, the skewness was less stable and less predictable. The observed abrupt surface warming on 16 July corresponds to stronger stratification, which is known to inhibit upwelling 2,4,24 . Therefore, the resulting vertical compression of the plume is inferred to have controlled this f-h Ortho-images of the site acquired by unmanned aerial vehicles (time is indicated by green crosses on the timeline): (f) before the deployment and calving event; (g) after deployment, but before the GLOF; (h) after deployment, during the GLOF (the colours are over-saturated to highlight the presence of the plume's core with a concentric turbidity pattern). i Helicopter view of surface ripples of the plume (not to be confused with calving or waterfall activity~70 m away, as known from direct observations and Supplementary Movies 1 and 3 response of the deep sensor. The exact mechanisms are difficult to constrain; however, we suggest that lowering of the plume top relative to the deep sensor (on a cable of a fixed length) somehow led to relatively stable gliding of the sensor in the plume.
It is known that ocean tides modulate ice motion and water circulation, however, there is a limited number of direct measurements near the calving fronts of marine-terminating glaciers in Greenland and Antarctica 21,22,[25][26][27] . Both the shallow and deep datasets were tide-modulated, as deduced from autocorrelation and wavelet-based semblance analysis comparing time-series with ocean tides (see Methods). The former was able to detect tidal modulation only in shallow-sensor data ( Supplementary Fig. 2), while the latter was particularly helpful for understanding deep-sensor data, allowing comparison of datasets of distinctly different types and accounting for the temporal variability of wavelength (Fig. 5). From direct GPS measurements, it is known that the terminus of Bowdoin Glacier is almost grounded 21,22 . For example, in July 2017, the terminus had a negligible tide-modulated vertical motion of <2 cm 26 (i.e. two orders of magnitude smaller than the tidal amplitude). Therefore, shallow signals were correlated with the tidal amplitude, owing to a deepening of the sensor into colder, saltier water at high tide as the water level rose around the sensor ( Fig. 5 and Supplementary Fig. 2). In contrast, the deep temperature data was negatively correlated with tidal rates of change, with temperature decreasing during rising tides, and vice versa. A weak positive semblance between tidal rate and pressure was also observed in the deep sensor data (Fig. 5), with the sensor probing deeper water during a rising tide, which we interpret as a dynamic deformation of the plume by the inflowing/outflowing tidal stream before and after a period of slack water.
Upward loops yield natural time-lapsed profiling of the water column. Ascents of the deep sensor eventually yielded timelapsed profiling of the water column, with 'fingerprints' of the plume's dynamics (Fig. 6). Specifically, we note (i) a downward shift of the cold-water layer (50-75 m depth) by~25 m during the period after the step-like surface warming (after 16 July), and (ii) short-term upward loops that are similar across different scales. The loops consistently exhibit recurrent features of ascent into cold, low-salinity water and descent into warm, saline water several minutes later. Such rapid and unexpected changes (1.0 ∘ C, 0.25 PSU) are due to entrainment of the sensor by a pulse of upwelling into Polar Water followed by the sensor's descent into water influenced by Atlantic Water 9,24 . Occasionally, the sensor moves via convection of water that is 0.5 ∘ C warmer than that of the same depths during descent (60-80 m), which explains the wide range in the temperature profiles (Fig. 3j). Sinking can be interrupted or slowed by entrainment into consecutive pulses. The recorded pressure oscillations and the observed inclination of the cable outward from the ice cliff characterised the motion of the deep sensor, with these being used to estimate current drag (see Methods). The equivalent horizontal flow rate necessary to produce the observed motion of the deep sensor is 1-2 m s −1 (Fig. 7).

State-space reconstruction by nonlinear time-series analysis.
From the preceding two sections we recognise that our next task is to reconstruct the high-dimensional nonlinear behaviour of the dynamic system from scalar measurements. A solution to this non-trivial problem of deciphering quasi-periodic and irregular Upper photographs show the lake before and after the GLOF; green rectangles depict the areas that were cropped and horizontally`stacked' for visual tracking of the water level (d). The seismic waveforms (c) and spectrogram (e) highlight the rise of low-frequency seismic energy due to intensified subglacial water discharge caused by the GLOF. Vertical red lines show the simultaneous onset and abrupt end of the increased seismic energy and water level change. Time is in seconds relative to 00:00:00, 13 July 2017 UTC. signals is a key question in nonlinear time-series analysis. An interesting feature of this growing discipline is 'time-delay embedding', which is used, for example, in electrocardiogram analysis (a brief introduction is provided in Methods). To unravel the complexity of data from both sensors, we unfold (or reconstruct) the time-series through time-delay embedding into the corresponding phase-state space. The essence of this framework is a graphical representation of the states of the underlying system (usually in 3D), which are best understood when animated (see Methods). We stress that this powerful tool requires a substantial effort on the part of the operator to visually inspect the dynamic features, and a major portion of what follows simply describes Supplementary Movies 5-9, which provide the most intuitive illustration of our analysis.
To provide a visual summary of the dynamics revealed by the five time-series (i.e. shallow-sensor temperature and salinity, and deep-sensor pressure, temperature and salinity), we were able to project our analysis onto a plane (Fig. 8) using the following steps: (1) we used a conventional representation of each time-series embedded in 2D, collapsing the 12 days of measurements together (upper subplots of Fig. 8a-e); (2) we created a 3D version of the same result, but using time as the z axis to reveal the extreme amount of otherwise hidden temporal variability in the behaviour of the system (lower subplots of Fig. 8f-j); and (3) to aid visual distinction of different water masses, we added the fourth dimension by using colour to indicate water temperature. This approach (detailed in Methods) allows visualisation of an overview of phase states, and of a particular state at a particular time across all five time-series. Hereinafter, by 'phases' or 'states' we refer to a set of conditions to which the system tends to converge, also called 'strange attractors' in chaos and nonlinear studies; and refer to notable jumps from one state to another as 'phase transitions'. Attractors are revealed by specific recurrent, but not exactly repeated, pathways of the system, which are also termed 'trajectories'.
The resulting shapes of the revealed phase-spaces are different among our time-series because each scalar is travelling through a unique scalar field, like pressure or temperature. For example, as detailed below, the deep pressure measurements trace a cone representing the sensor's motion (Fig. 8a); in contrast, the deep temperature data trace a square due to travelling of the sensor between two different water masses with varying duration of stay within each of them (Fig. 8b).
The phase-space trajectories of the deep sensor data. Even when the deep sensor is episodically disturbed, the pressure values exhibit a cone-shaped region of the phase space, where close-by trajectories (i.e. recurrent orbits) converge (Fig. 8a). The tendency towards this specific state (i.e. a strange attractor) is due to a balance between the downward gravitational force and the upward drag of the plume and icebergs, which occasionally tugged on the sensor's cable. The state-space trajectories are smooth, as expected for sufficiently sampled sensor motion underwater for the maximum number of significant changes set to three). Dashed horizontal lines show means and associated differences, ΔT, for the segments before/after the afternoon of 16 July 2017 for water temperature, air temperature, and minimum night air temperature (blue circles).
(Supplementary Movie 5). In the evening of 14 July, the onset of perturbations leading to temporary divergence from the attractor coincided with the GLOF event ( Fig. 8f; Supplementary Movie 5). At the end of the GLOF, the trajectories settled into a deep, calm vertex of the attractor for approximately half of a day (Fig. 8a), before being quickly ejected up to the water surface followed by a continuation of wide-loop, self-similar trajectories ( Fig. 8f; Supplementary Movie 5). The latter prolonged ejection, with a square-like trajectory, was most likely due to an iceberg pulling the cable (14:38-16:17, 15 July). The amplitude of perturbations gradually decreased towards the end of the observations (Fig. 8f).
Reconstruction of a state space for deep-sensor temperature (Figs. 8b and g; Supplementary Movie 6) showed that the initial state (~−0.6 ∘ C) was perturbed by the GLOF, which pushed the sensor upwards into slightly colder water (~−1.0 ∘ C). After a transitional period of~9 h, the sensor reached deep water with the highest temperature for this depth (~−0.4 ∘ C; Figs. 8b and f), before a change of state back to the initial attractor, but from then on displaying occasional loop-like (10-15 min long) fluctuations into colder (but shallower and less dense) water (~−1.5 ∘ C).
Time-delay embedding of deep-sensor salinity showed the most drastic phase transition due to the GLOF (Figs. 8c and h; Supplementary Movie 7). First, the state of the initial small attractor (33.7 PSU) was destabilised. Second, after the end of the GLOF, after an elapsed time of 9 h, the system transitioned through a narrow bottleneck into the second attractor of lower salinity (32.7 PSU; termed 'post-GLOF anomaly' and marked with a circle in Fig. 8c), corresponding to the greatest depth (Fig.  8a). Afterwards, the initial state was then re-established, but was frequently perturbed by brief departures of~−0.2 PSU from the centre, characterised by oscillations between a warm, saline attractor and a cold, less-saline repeller; i.e. an anti-attractor or set of conditions where the system does not 'want' to remain (Fig. 8h). During the prolonged phase transition and orbit within the secondary attractor (Fig. 8c), the sensor probed fresh water during occasional ascents (limited to depths of >80 m), at depths that had higher salinity before and after this anomaly (except brief passages of similar water during manual profiling; Fig. 3k and Methods). This suggests that after the GLOF, the sensor was most likely trapped in an ice-cliff cavity between the ice and the upwelling core with an anomalously high ratio of freshwater (~3%) for such depths 24 , until it was released by the aforementioned iceberg. Compared with pressure and temperature, the temporal evolution of salinity exhibited trajectories that were less smooth, although even with jump-like dynamics the data points recorded within a 1-minute interval are clustered together (Figs. 8c and h; Supplementary Movie 7).
The phase-space trajectories of the shallow sensor data. The embedded shallow salinity and temperature sensor showed less evidence for perturbation by the GLOF and larger amplitudes of variability relatively to the deep sensor (Figs. 8d, e, i and j); however, a slight decrease in temperature and concurrent increase in salinity were observed during the GLOF (Supplementary Movies 8 and 9). Similar changes could be produced by a rising tide by deepening the sensor into slightly colder and saltier water. However, the observed changes took place during the falling tide, which should produce the opposite effect (Fig. 5), suggesting that the GLOF brought entrained cold and saline water to the surface, as previously observed at Narsap Sermia Fjord 13 .
The most significant departure from the initial state is observed in the afternoon of 16 July, when temperature rose from about 0 to +2 ∘ C (Figs. 8d and e) and salinity dropped by 3 PSU (Figs. 8i and j). We note that after this departure, the deep sensor maintained a negative skew in the rate of change in pressure ( Supplementary Fig. 1c). More extreme intrusions of fresh water began after 20 July, and early in the morning on July 23 the system moved into a new highly unstable, low-salinity attractor (25 PSU, +5 ∘ C; Figs. 8d and e). The presence of such warm water at this shallow depth indicates that after 17 July, near-surface fjord conditions dominated the behaviour (Figs. 8i and j). The only exception was observed after an episode of strong northwest wind (18 July, 14:30-18:30; Supplementary Movies 3-4 and Supplementary  Fig. 3). A few hours later, at around 22:00, ascent of the deep sensor up to almost 30 m depth was followed by a 2 ∘ C decrease in temperature and a 1.5 PSU increase in the shallow-sensor salinity lasting at least 15 h (Figs. 8i and j). This transition could be seen in the time-series (marked with circles in Figs. 3g and h), and led to an increase in density (Fig. 3i). It was also associated with a change in skewness and kurtosis of the shallow salinity observations (i.e. statistical measures describing the shape of a distribution; Supplementary Fig. 1c). This suggests a temporarily weakened stratification on July 19, which was likely preconditioned by a wind shear stress at the water surface (for more details on the potential wind forcing, refer to a phenomenon of upwelling due to the Ekman drift 28 ).
Power spectra and energy cascade. Time-series spectra were highly variable (Fig. 9). All deep sensor spectra were relatively smooth, whereas the shallow sensor spectra displayed a broad peak around 10 min with a steeper roll-off at high frequencies.
For deep signals, a strong rise and steepening of spectra could be observed for all tracers during the GLOF. Deep pressure spectra were the steepest and lacked energy at periods shorter than~7 min, presumably because of the dampened response of the system underwater (e.g. drag on the cable might have limited the ability of the sensor to respond promptly).
Comparison of the global median spectra with particular episodes (i.e. before/during/after the GLOF), demonstrates the difficulty in understanding the activity of the plume from spectra alone. This is particularly true for the shallow signals, which were dominated by surface processes caused by lateral transport of water along the front, when the sensor might not be reached by the plume. For example, we sometimes observed a counterclockwise gyre in front of the plume (15 July; Supplementary Movie 1), indicating that horizontal-shear-induced turbulence could influence shallow spectra.
The presence of a power-law scaling and approach to the -5/ 3 spectral slope, especially for the deep-sensor instantaneous conductivity measurements (Fig. 9), is in line with the Kolmogorov turbulence theory for a so-called inertial subrange of a direct energy cascade 29 . The cascade represents scale-free transmission of externally injected energy between large-scale structures (or eddies) and smaller structures, and is bounded by a so-called dissipation range, where mechanical energy is converted to thermal through viscous molecular friction. The shallow signals, however, could be influenced by mechanisms such as shear near the cliff and 2D turbulence, which is known to form a split cascade with energy transfer in both directions 30 . A lack of intermittency and a Gaussian distribution are properties of an inverse energy cascade in 2D turbulence 30 , which we, however, did not observe ( Supplementary Fig. 1).
Autocorrelation as a measure of the 'memory' of the flow. Since the passage of disturbance (e.g. a cyclone) over a receiver takes time, signal autocorrelation can be considered a measure of the 'memory' of the flow (see Methods for details). The dynamic autocorrelation structure of the data shows correlation up to 300 s (Fig. 10) and highlights the occurrence of coherent events. Median integral time scales,τ (defined by Methods' Eq. 2), were between~28 and 140 s. For the deep sensor, the longestτ were observed for pressure measurements (140 s), the shortestτ were typically observed for salinity (28 s), and the temperature was 110 s. For the shallow sensor, in general, variability in integral timescales was narrower than that for the deep sensor. Moreover, in the shallow sensor, theτ of salinity (48 s) was also shorter than that for temperature (66 s). Integral-time-scale separation between slow temperature and fast salinity is pronounced at the deep sensor, whereas the shallow data do not exhibit such separation.
The difference between shallow versus deep autocorrelations suggests that individual flow structures could not be detected upstream because they either evolve or did not travel between the two sensors (e.g. the halocline at~15 m depth could have been a barrier).
In general, using a comparison of the statistical properties of the dominantτ, pressure appears to be a slowly changing process while salinity appears to rapidly change. This result is reasonable because the motion of the sensor in the fluid is smooth due to drag forces, whereas turbulent salinity microstructures (with sharp gradients) can pass through the instrument much more quickly. Owing to the high-frequency variability (also seen in our PSD results as shallower spectral slopes; Fig. 9) the system 'forgets' its previous state quickly.
According to Taylor's hypothesis of frozen turbulence 31 , the largest integral timescales can be related to length scales using the average flow velocity, u, as L ¼τu. This assumption about flow advecting the turbulence is the only possibility for estimating spatial information from single-point time-series data and may not be valid in all cases 32 . For a sensor dragged through a passing fluid with uncertain velocity within a layered structure, the above assumption is not trivial, as it may correspond to modified correlation lengths. Taking into account these caveats, and assuming a velocity of 1.5 m s −1 , the dominant length scales of coherent disturbances are between~40 and 210 m, depending on the scalar. The choice of velocity is necessarily arbitrary, although, it is generally representative of inverted current speeds (Fig. 7), and modelling studies have used values of 0.3-3.1 m s −1 at the source of a plume 16 . However, the estimated length scales are reasonable as the largest eddy scale is expected to be limited by the water depth (~200 m).
Furthermore, autocorrelation results show the existence of oscillatory features (Fig. 10): where the shape of the autocorrelation after the first zero crossing may have temporally persistent structures repeating periodically. This could represent a population of similar-scale eddies passing through the sensor in a sequence train. Such eddies play a dominant role in the entrainment process by engulfing external fluid 33 . Implications, conclusions and outlook. This study provides the first example of continuous oceanographic monitoring in the most difficult of cryospheric environments. Calving fronts remain among the least-explored oceanic zones, with the main implications of the study being instrumental, methodological, and physical, as follows.

To obtain in situ observations at the calving front of a
Greenlandic marine-terminating glacier, we pioneered a novel instrumental approach and incorporated our observations into the most comprehensive plume-monitoring campaign to date (Fig. 1b), including the first simultaneous monitoring of both water intake and output (lake and plume, respectively). This high-temporal-resolution dataset may become an important reference for numerical simulations and follow-up studies. 2. We developed new methodology for deciphering oceanographic Lagrangian records of convoluted fluid and solid interaction through nonlinear time-series analysis (i.e. time-delay embedding), allowing observation and characterisation of a deterministically chaotic system (Fig. 8). This framework is unconventional in fluid mechanics, geophysics, and oceanography, but has the potential to become a powerful analytical tool in geoscience. 3. Our sensor data provide a complex portrait of plume and fjord dynamics with a rich diversity of processes operating on short time-scales, which have been undocumented to date (as detailed below). For example, a sudden eruption of a subglacial discharge plume is a low-frequency co-seismic process (Fig. 2) accompanied by intensified water-surface ripples (see Supplementary Movie 1 or Fig. 1i illustrating the general view of such ripples), with seismic and seasurface-height observations thus opening new avenues for better understanding of the plumes.
To contextualise our observations for submarine melt, we note that previous studies of ice-ocean interaction have been based on two fundamental elements: (1) limited and rare conductivitytemperature-depth (CTD) casts representing background stratification in models; and (2) subglacial discharge calculated from runoff simulations 4,5,16,24 . These two fundamental elements have been indispensable for learning a lot about the first-order dynamics of ice-ocean interactions, especially in the long-term. Our study results indicate that, at least in the short-term, neither of these elements is completely valid, because stratification can change within a few days. A single CTD cast per summer season provides a limited representation of conditions, and glacial-lake outburst floods can cause major disruptions to stratification. This is of importance to our understanding of submarine melting. For example, even if a single CTD cast could provide a reasonable representation of the slowly changing deeper water masses throughout the time period, an abrupt~25-m-downward shift of stratification (Fig. 6a) means that, near the surface, the contact area between the relatively warm water of the fjord and the ice cliff increases by~70,000 m 2 along the 2.8-km-wide terminus. It also raises the question of how such a sudden shift enhances  Spectral analysis. Power spectral density computed by using 2-h-long sliding windows band-pass filtered between 12 h and 20 s. Thick black curves indicate the median spectra of each time-series (also assembled in the lower-right plot where solid and dashed curves correspond to the deep and shallow sensors, respectively). Thin black lines show spectra before the GLOF (14 July 00:00-02:00), during the GLOF (July 14, 23:00-July 15, 01:00), and after the GLOF (July 19, 00:00-02:00) during a pulse of cooling and salinity increase observed near the surface (the legend is provided with`Shallow Salinity' subplot). Diagonal lines show the Kolmogorov inertial subrange spectral slope of −5/3 for reference. The`Nyquist' and`Aliasing' blue vertical lines are related to the sampling interval (Δt) and show the frequencies below which the spectra are correctly resolved (the Nyquist frequency, 1/(2Δt), and thè Aliasing' frequency, 1/(4Δt), respectively).
undercutting (i.e. ice-cliff overhanging) and impacts stress distribution associated with crevassing and calving. This highlights the necessity of further continuous subsurface observations for verifying how wide-spread and frequent such changes are, and what their long-term impact is. Other questions arise from evidence of the tidal modulation ( Fig. 5 and Supplementary Fig. 2). For example, the influence of tides on plumes and melting does not seem to have been addressed in Greenland, but is recognised as an important driver of subaqueous melt for ice shelves fringing Antarctica 27 . Tideinduced currents and vertical motion of background stratification may deform the plume and modulate local melt rates, requiring further in-depth studies.
In addition, we suggest that the plume is also influenced by wind forcing (Figs. 3g-i and Supplementary Fig. 3). This further highlights that current models of convective plumes are lacking some processes. Considering that the ocean is strongly influenced by wind shear near the surface, the role of wind remains to be understood.
At fine scales, our spectral-analysis results (i.e. tracer variance) are consistent with a direct cascade of energy (Fig. 9). Turbulence has long fascinated scientists searching for its universal properties, one of which is the cascade. We acknowledge that, because of the buoyancy and shearing of rising water, high-energy turbulence in the plume was expected and assumed by the research community. However, there has been little in situ evidence that characterised the plume as a turbulent flow. Ocean turbulence is generally fundamental to heat and matter transport and has different regimes, many of which remain unexplored and are difficult to assess. Our results imply that the celebrated Kolmogorov-1941 framework 29 extends to turbulent plumes, which is a far removal from the idealised turbulent flows for which the framework was originally postulated.
On the one hand, the spectral statistical similarity to conventional turbulent flows is an aid to parametrising subglacial discharge plumes in numerical models; on the other hand, in addition to GLOF-like runoff events, the occasional current pulse means that the assumption of subglacial discharge being a steadystate process might not hold, at least during the period of our observations.
The lack of direct velocity measurements means it is difficult to use our experimental results to evaluate previous models, particularly those involving turbulent entrainment parameterisations. Buoyant-plume theory is intended mainly to represent the broad outlines of a convective plume dynamics without resolving the mechanics of individual turbulent eddies being necessary 1 . Such turbulent gravitational convection theories are relevant to plumes that have reached a steady state. The high-energy, intermittent turbulence shown here suggests that such a state could be unattainable. We also note that time-lapse photography of the plume area (Supplementary Movie 1) shows an overwhelmingly complex high-frequency temporal variation in conditions near the calving front. This implies that some dynamical effects, such as those caused by episodes of intensified entrainment, might be overlooked. It may, therefore, be unrealistic to assume a certain entrainment rate, with only numerical experiments involving an intermittent buoyancy flux from the source clarifying impacts on a plume's long-term role in a fjord. There is also a broad spectrum of processes missed by simplified numerical treatment of the plume. For example, turbulent structures are known to continuously shape the geometry of rivers 34 , and this should also apply to ice-cliff morphology. The efficiency of the plume in melting the ice (icecliff face, sea ice, icebergs, and ice mélange), entraining water, and transporting suspended matter to the surface strongly depends on turbulence. Furthermore, a growing body of literature highlights the importance of turbulence to marine organisms and ecosystems. For example, turbulence induces active adaptation of zooplankton to control their distributions and the dispersal of their populations 35 , and evolution of chemotactic algorithms for plume-source tracking by many animals 36 . Only eddy-resolving models can be used to study such processes in detail.
In the light of the evidence presented here and discussion of the potential impacts of the high-frequency variability on submarine melting and material fluxes, the data demonstrate a need for the consideration of intermittent turbulence, day-scale variability in stratification, and tidal/wind effects in models developed to quantify nutrient fluxes and plume-driven melting, which are of major importance to biogeochemistry, glacier retreat, and sea-level rise 4,37 . We acknowledge that there is a need to understand how (or if) such high-frequency variability influences the longer-term mean processes and how this behaviour could be realistically included into large-scale models of the ocean or ice sheet.
We anticipate our observations will constitute an urgently needed first step towards establishing constraints for nextgeneration eddy-resolving models, which are required for understanding key processes such as glacier sub-aqueous melt and water circulation and mixing, as well as the effects on marine biology. Hopefully, such work will increase our ability to reliably predict the response of tidewater glaciers to warming, estimate their contribution to sea-level rise, and forecast the evolution of fjord ecosystems.

Methods
Seismic and lake measurements. To record glacio-hydraulic tremor, we installed a seismic station directly on the ice, about 300 m up-glacier from the plume. A three-component Lennartz LE-3D short-period sensor (eigen-period of 1 s; flat response between 1 and 100 Hz), connected to a DATA-CUBE 3 recorder, a battery, and a solar panel, was placed in a shallow ice pit on a metal tripod and covered with a high-albedo protective blanket. The station operated with a sampling frequency of 400 Hz. The instrument-corrected signal was converted to velocity and then bandpass-filtered (0.3-1.0 Hz) with a zero-phase-shift, 4th-order Butterworth filter. The spectrogram was constructed using 1-h-long windows with 50% overlap.
We recorded temperature and water pressure every minute using a sensor (HOBO U20, Onset Co.) installed at the coast of the ice-dammed lake approximately 420 cm below the lake surface. An impulse wave generated by calving of marginal dead ice into the lake on 14 July (00:41) displaced the sensor. Near the end of the drainage of the lake (14 July, 23:25), receding water stranded the sensor. We converted the data to water level by removing the standard atmospheric pressure, P SAP and converting the remaining pressure to equivalent water height (using h w = (P Total − P SAP )/ρ w g, where ρ w is the density of fresh water and g is the acceleration due to gravity).
Time-lapse photographs (4000 × 3000 pixels) of the ice-dammed lake were taken using a GARMIN VIRB ® XE camera with a 2.62 mm lens and GPS clock. The camera was installed on rock, powered by a battery and a solar panel, and captured images every minute. Vertical slices of an ice cliff (60 × 100 pixels) were cropped from 4,320 photographs (taken between 00:00, 13 July and 00:00, 16 July) and aligned to complement the missing pressure-sensor data (from sensor displacement discussed above) and highlight the overall water-level change associated with drainage.
Fjord observations. Conductivity plus temperature and pressure were measured with a conductivity sensor (4419RB) and a pressure sensor (4117BR), respectively (by Aannderaa Data Instruments, Norway). Both sensors had RS-422 full-duplex interfaces for use with long cables. Sensors were anchored by~5 kg weights (rockfilled bags) and connected with cables (diameter: 8.3 mm, weight: 0.1 kg m −1 ) to the SmartGuard recorder (Extended version 5300 by Aannderaa Data Instruments), which was powered by an external battery and a solar panel. Resolutions of conductivity, temperature, and pressure were ±0.0002 S/m, ±0.01 ∘ C and ±0.0001% of the full scale output (FSO; 0-4000 kPa range), respectively. Accuracies of conductivity, temperature and pressure were ±0.0018 S/m, ±0.05 ∘ C and ±0.02% of FSO, respectively.
The pressure sensor had an instantaneous response time. The response time (90%) of conductivity measurements was less than 3 s, whereas the response time (63.2%) of temperature measurements was less than 10 s. In other words, relative to the adopted sampling rate of 10 s, pressure and conductivity can be seen as fastresponse, instantaneous sensors; the thermistor, however, could have a lag of a few tens of s and therefore a lack of energy at frequencies higher than 0.02 Hz due to the instrument response. This mismatch in the response of conductivity and temperature measurements implies that the calculated salinity inevitably becomes their composite 38 . However, we note that heat has higher diffusivity than salt by two orders of magnitude (different Batchelor scales) 39 , which may correspond to smoothing of temperature gradients but sharper salinity features, and thus may lead to fast, non-instrumental disappearance of high-frequency content in the temperature spectra.
On 24 July 2017, a conductivity-temperature-depth (CTD) cast was collected using an ASTD 102 JFE Advantech profiler from a boat in the middle of the fjord, 1.4 km from the terminus ( Fig. 1b; see 9,40 for sampling details).
Two unmanned aerial vehicles (UAVs) were used ( Fig. 1f-h). The vertical takeoff and landing, fixed-wing UAV 'WingtraOne' (www.wingtra.com) was used to conduct a photogrammetric survey on 7 July. The fixed-wing UAV 'SenseFly eBee' (www.sensefly.com) was used to conduct surveys on 13 and 14 July. The UAVs were both programmed to take aerial photographs with an overlap of 90% in flight direction and 75% in cross-flight direction. Ortho-images of the calving front with a resolution of 0.5 m were processed from UAV aerial images by structure-frommotion photogrammetry using Agisoft PhotoScan software (www.agisoft.com).
Time-lapse photographs (7952 × 5304 pixels) of the calving front were taken every minute from Sentinel Nunatak (Fig. 1) by a Sony ILCE-7RM2 camera with a Sony FE 28 mm F2.0 lens and an external GPS clock (4-17 July). To zoom into the location of the plume set-up, sections of the photograph [1100 × 367 pixels] were cropped. Additional time-lapse photographs (1280 × 1024 pixels) of the plume area were taken every 5 or 15 min from the southeast side of the fjord by two co-located GardenWatch cameras (Brinno Co.) with 5.01 mm focal length and a hand-watchset clock (time precision ±30 s; 16 July to 01 August). Photographs were synchronised with deep-sensor pressure, shallow-sensor temperature, tide, and air temperature, and were compiled as Supplementary Movies 1, 3, and 4. To highlight changes, each image was converted to grey scale and subtracted from the previous image.
Conversion of data. Strong tidal modulation of the shallow sensor signals (see below) suggested variable depths of measurement caused by water-level oscillation above the sensor. To estimate the variation of depth with time, we modulated the initial deployment depth (5 m) by the tidal amplitude measured at Thule (the same phase and amplitude as at Bowdoin Fjord; 22,41 ), yielding depths between 4 and 7 m.
Eulerian representations of the Lagrangian dataset are provided in Fig. 3g-h and j-k. In addition, we calculate water salinity and density 42,43 . Density is controlled mainly by salinity and density variability corresponds to time-varying stratification and noise (Fig. 3l). In particular, we note that for the deep sensor, the cross-correlation identified a 10-20 s (1-2 samples) lag of temperature behind conductivity owing to a slower thermistor response at cold temperatures and high pressures ( Supplementary Fig. 4). Non-consideration of this lag leads to a mismatch and numerous, spurious spikes in default-recorder-computed salinity exceeding 34.2 PSU (pressure set to 0 kPa), which could be misidentified as entrained Atlantic Water. To address this issue, we adopted a shift in temperature by -20 s as it reduces high-frequency pollution in computed salinity, and we used this shift in computation throughout the paper. Furthermore, we consider timevarying pressure in all salinity computations.
Autocorrelation. The autocorrelation coefficient corresponds to a normalised correlation of a variable with itself and is defined by: where τ is the time lag, σ 2 is the variance (the root-mean-square of x), and 〈 ⋅ 〉 represents averaging. According to Schwartz' inequality, |R| ≤ 1 for all τ.
If the sampling rate is sufficient, eddies passing through a point correspond to a non-instantaneous change, and thus the signal remains correlated to the previous time-step. For example, large structures, like synoptic-scale eddies, need a long time to pass by 44 . Therefore, the rate of decorrelation can be seen as a measure of the 'memory' of the flow. Fluctuations in random (or under-sampled) signals decorrelate quickly with the lag τ.
The integral time scale of a variable is: We follow 45 and take T corresponding to the first zero crossing (i.e. R (t, T) = +0).
For a short-term autocorrelation, we use a demeaned data segment 32 within a 30-min-long sliding time window with a highest possible increment step of 10 s (i.e. leading to an overlap of almost 99.5%); the maximum lag is limited to ±30 min. For a long-term autocorrelation (as a way to identify tidal signal), we filter out periods shorter than 100 min and use a window 48 h long with a 10 min increment (99.6% overlap); the maximum lag is limited to ±48 h.
Continuous Wavelet Transform. The continuous wavelet transform (CWT) of time-series data allows analysis of temporal changes in the frequency content with higher resolution compared with the Fourier transform. The CWT of a time-series, x(t), is defined as 46 where u is displacement, s is scale, Ψ is the mother wavelet used, and * corresponds to a complex conjugate. Following 46,47 , here, the time-series is convoluted with a scaled complex Morlet wavelet, which is given by where f c is the wavelet centre frequency, and f b is controlling the wavelet bandwidth. To compare the continuous wavelet transform of the ocean tide amplitude, T, or its derivative (CWT T or CWT dT/dt ) with our time-series (CWT i ), we first resample and interpolate the tidal data (5 min sampling rate) to a common time frame with the rest of the data, and then multiply an amplitude of the complex cross-wavelet transform (CWT T,i = CWT T × CWT Ã i ) by a function of a local phase given by where the cos(. . . ) part measures phase correlation between two inputs and represents the wavelet semblance, and D represents a dot product de-noising the latter 46 . Through this nonorthogonal wavelet analysis, we find that the tidal signal in the 'deep' data correlates with the rate of water-level change (dT/dt), whereas that of the 'shallow' data correlates directly with the tidal amplitude (Fig. 5).
Inversion of current velocity. Variation in pressure (Fig. 3d) and inclination of the cable outward from the ice cliff observed from a helicopter (17:56, 17 July; Fig. 7) indicate a pendulum-like swing of the deep sensor on a line of a fixed length, L, due to the strong outward push of the plume. This configuration is similar to a reverse problem of current drag and tilt of buoyant oceanographic moorings 48 . To estimate the horizontal current speed needed to produce the observed pressure variation, P(t), we make the following assumptions.
First, we neglect any drag on the thin line as a secondary effect and focus on the drag force, D f , on a ballast at the end of the line, given as 48 : where ρ w is the average density of water (1,027.3 kg m −3 ), C D is the drag coefficient, A is the frontal area of the body, and v is the velocity of the current flow, which we want to calculate.
Second, we assume that to sustain a certain angle, α, the weight has to be pushed against the gravitational pull. The angle is where I c is the ice-cliff height (30 m), and the length of L = [max(P(t)) + I c ] (corresponding to the deepest and 'quietest' conditions on 15 July).
The effective gravity underwater is where ρ r is the density or rocks (2,700 kg m −3 ), and g is the acceleration due to gravity (9.8 m s −2 ).
For a certain instantaneous α(t i ), the equivalence of the tangential force components of current push, F ct , and gravity pull, F gt , implies that the horizontal pushing force, F c , is Equating F c and D f yields To address poorly constrained parameters (m, C d , A), we compute histograms of v x (292 bins) for an ensemble of 224 parameter combinations, where m ∈ [4, 10] kg, C d ∈ [0.5, 1.2]ref. 48 and A ∈ [0.0157, 0.0289] m 2 and [0.0595, 0.0893] m 2 (the former area is consistent with the cross-section of a sphere of a given mass, whereas the latter is the increased area approximately corresponding to the folded bag used for rocks).
Embedded state-space trajectory. A power spectrum is insufficient either for describing the complexity of nonlinear time series, as two signals with identical spectra may have drastically different dynamics 49 ; or for detecting the presence of attractors 50 . To overcome this, so-called 'time-delay embedding' was proposed for state-space reconstruction of nonlinear dynamical systems using one observable 50 . Mathematically, this revolutionary nonlinear-time-series analysis approach is a diffeomorphism, based on a system of delay coordinates, which allows unfolding (or reconstructing) the attractor that is topologically similar to the true attractor of the original dynamical system 50 .
To properly reconstruct the attractor from a scalar time series, it is necessary to choose two key parameters, namely, the time delay, τ, and the embedding dimension, m e . Currently, the optimal determination of these two parameters remains mathematically non-rigorous 49 . The principal requirement is to choose such a time delay that the consecutive coordinate is independent. In linear series, such delay corresponds to the decorrelation time (i.e. the first drop of the autocorrelation function either below e −1 or zero), whereas in nonlinear analysis it is common to rely on the self-mutual information criterion I(T), which corresponds to the amount of information we learn about x n by measuring x n+T 49,51 .
The embedding dimension is typically determined by false nearest neighbour (FNN) analysis. This technique differentiates between genuine and 'false' neighbours, where genuine neighbours remain close in the higher-dimensional embedding (R m+1 ), and the proportion of false neighbours (which could be geometrically close simply because of insufficient unfolding in R m ) goes to zero. In practice, m e is the smallest value at which the fraction of FNN(m) lies below some small arbitrary threshold (or zero) and is typically fixed as 3 or 4 49,51 .
I(T) can be sensitive to histogram binning (b) while estimating a probability distribution, whereas FNN(m) is a function of the chosen τ and varies with the socalled ratio factor (R), which depends on the spatial distribution of the embedded points 51 . Following sensitivity tests, we found that b was of little importance to the shape of the I(T) curve; b was fixed to 2N 1 3 , where N is the number of points, and R was set to 20, as commonly used 51 .
Both I(T) and FNN(m) were computed for 12-h-long segments of all data and are shown in Supplementary Fig. 5. The results suggest that for most of the data, τ plateaus above~5 min for most observables, for which the data are appropriately 'unfolded' at m e = 3 − 4. Deep salinity shows the lowest self-mutual information, which flattens above~2 min ( Supplementary Fig. 5); similarly, it has the shortest decorrelation time ( Fig. 10; in particular, deep-sensor salinity fluctuations are the fastest process, and occasionally have a near-zero integral timescale). We suggest that the presence of the aforementioned noise in the deep salinity data leads to deviations in both I(T) and FNN(m). Nevertheless, sensitivity embedding tests with the lower embedding dimension of 2 and shorter/longer delays (~2-16 min) yielded topologically equivalent results. The smaller dimension implies that instead of the attractor, the reconstruction represents its projection. Therefore, for the sake of simplicity and consistency, here we present 2D embedding with τ = 5 min (Fig. 8). Animated results giving intuitive understanding of the system state dynamics are provided in Supplementary Movies 5-9. To fully reveal the temporal evolution and state-space architecture of the system, we plot the 2D embedded time-series ( Fig. 8f-j) as a function of time (Z-axis) and the corresponding temperature (by colour).
Spectral analysis. In general, turbulent structure and statistics are inferred from Reynolds decomposition of the velocity signal 52 . Such data can be used either for estimating turbulent kinetic energy dissipation rates from power spectra (i.e. among the most needed for plume models; 2 ), or for finding the characteristic size of the largest energy-injecting eddies from autocorrelation 36 . In our logistically challenging case, velocity could not be measured directly. As an alternative to in situ measurements of the turbulent velocity field [e.g.39,52 , quite commonly in experimental studies, microscopic tracer particles or hydrogen bubbles are added for analysis of turbulent flow via particle image or laser doppler velocimetry 36,53 . Acoustic scattering techniques quantify turbulent ocean environment in which either biological or physical sources of scattering are present 54 . Other passive and active tracers carried by the flow, such as temperature and salt, can be used as a proxy for turbulence. This has been recognised by theoretical work showing that tracer variance has spectra similar to that of the flow 36,55-57 . Tracer variance has been successfully applied in oceanographic studies, including an investigation that found that salinity variance may sufficiently characterise turbulent fluctuations for cases where salinity dominates density, as in our case 39 .
Therefore, we compute the one-sided power spectral density (PSD) of the timeseries that were bandpass-filtered between 12 hours and 20 s (Nyquist frequency) using fast Fourier transform and 2-h-long sliding time windows overlapping by 50%. Spectra were smoothed with a coefficient of 10 after 58 . For separating the possible influence of temperature-sensor inertia on salinity, in spectral analysis we also consider conductivity separately.

Data availability
All oceanographic, meteorological, and lake-pressure measurements are provided with the paper (Supplementary Data 1, Supplementary Movies 2, 5-9); these data are also publicly available at https://zenodo.org/record/4552842 59 . Freely available SENTINEL-2A satellite imagery was downloaded from https://earthexplorer.usgs.gov/. Thule tide data are publicly available through the Global Sea Level Observing System network (http://www.ioc-sealevelmonitoring.org/map.php). Seismic and time-lapse data for the lake are publicly available through the Arctic Data archive System website (https://ads. nipr.ac.jp/dataset/A20200108-005). Meteorological data are also described and partly available at https://ads.nipr.ac.jp/dataset/A20181016-004. Other time-lapse and UAV imagery is presented in figures and Supplementary Movies 1, 3, and 4. All Supplementary Movies of the paper (1-9) are also provided in high resolution at https://doi.org/10.5281/ zenodo.4552847 60 .