Detecting Thermohaline Layers with a Clustering Algorithm (under construction)

A few hundred meters below the surface of the Arctic Ocean, there is a layer of Atlantic water which is warmer than the top layer in contact with the sea ice above [10]. If all the heat in this Atlantic layer were to rise to the surface, all Arctic sea ice would melt within five years [11]. Winds across the ocean's surface create downward-moving internal waves which can cause vertical mixing of heat [3]. However, in the top of the Atlantic layer, the water's density changes suddenly every few meters forming a density “staircase,” a barrier to vertical mixing at each step [8]. As Arctic sea ice declines, it will be important to predict how the interaction between internal waves and density staircases might change.

Working with Prof. Nicolas Grisouard, I used numerical simulations to investigate this phenomenon and successfully reproduced the results of a wave tank experiment done in Prof. Peacock’s ENDLab at MIT [4]. These simulations give insight into whether the propagation of internal waves through density staircases could reasonably become significant to vertical heat transport as the Arctic climate continues to change. This work is presented in full detail in Chapter 4 of my thesis.

This page gives a brief overview of my work with my collaborators, Prof. Erica Rosenblum, Dr. Jonathan M. Lilly, and Prof. Nicolas Grisouard, on using a clustering algorithm to identify thermohaline layers in hydrographic profiles from the Arctic Ocean.

Background

A couple hundred meters below the surface of the Arctic Ocean, there is a layer of relatively warm water originating from the Atlantic Ocean [10]. This layer is evident in Figure 1, which shows typical temperature and salinity profiles measured by Ice Tethered Profilers (ITP’s) in the Canadian Basin [7], as the temperature maximum around 350 meters depth. In the ocean, 1 decibar (dbar) in pressure is equivalent to 1 meter in depth to a very good approximation [x]. ITP’s have a vertical resolution of around 25 cm [10].

The Atlantic layer contains enough heat to melt all Arctic sea ice if it were somehow able to rise to the surface [11]. However, the Arctic Ocean is primarily stratified by salinity as opposed to most other oceans which are stratified by temperature [x]. Therefore, because the Atlantic layer is relatively salty, it sits stably beneath the relatively fresh surface waters which generally come from river runoff and ice melting [x].

Just above the Atlantic layer, there generally exist stratification structures called thermohaline staircases, an example of which can be seen in the insets of Figure 1.

Figure 1: The temperature profile number 1257 from ITP 1 (left) and the corresponding salinity profile (right) taken on Monday, June 26th, 2006 at 136.1525◦W, 77.2870◦N. Both profiles show insets of the same depth ranges where double-diffusive staircases can be seen. Following Shibley et al. (2017) Figure 3b [7].

Figure 2: Reproduced from Schmitt et al. 1987, Figure 3a [x].

Thermohaline Staircases

Many areas of the World Ocean contain thermohaline staircase structures, also sometimes known as double-diffusive staircases [vdB]. These consist of a series of layers, where temperature and salinity are vertically well-mixed, separated by interfaces where the temperature and salinity change rapidly. The layers are generally 1—2 orders of magnitude thicker than the interfaces.

These structures have been shown to stretch horizontally for hundreds of kilometers [x]. This gives the stair-steps aspect ratios of 10⁵—10⁶ and are therefore sometimes referred to as “sheets” [x]. Figure 2 shows a series of profiles taken along a ship track and shows how, given the profiles are taken within a certain distance of one another, you can visually trace certain stair-steps across different profiles covering many kilometers in distance.

Thermohaline staircases were first observed in the Arctic Ocean since 1969 [x]. Since then, they have consistently been present in every major observation campaign in the Arctic Ocean including AIWEX and SHEBA [x]. While such campaigns have shown that thermohaline staircases structures are consistent features, to the best of my knowledge, no study has yet to show whether individual layers of staircases remain coherent over periods longer than a few years [x].

Detecting thermohaline staircases in data

too much data to look manually

ITPs have taken many profiles

other studies have come up with algorithms to identify staircases, but all work on individual profiles (show image from van der Boog’s paper which shows the selection of layer points and interface points)

I was inspired by timmermans 2008 figure below to try and use a clustering algorithm to automatically connect points within the same cluster across profiles.

Figure 4: The values of potential temperature and salinity for all points within the up-going profiles taken by ITP2. Reproduced from Timmermans et al. 2008, Figure 5a.

The HDBSCAN Clustering algorithm

Lots of other examples of clustering algorithms being used in oceanography.

I choose to use HDBSCAN for a number of reasons.

I will present the basic principles here, but would suggest looking at:

the original paper Campello, R., D. Moulavi, and J. Sander (2013) “Density-Based Clustering Based on Hierarchical Density Estimates,” Lecture Notes in Computer Science, doi: 10.1007/978-3-642-37456-2_14

the python package’s docs page: How HDBSCAN Works

these articles from towards data science:

“Understanding HDBSCAN and Density-Based Clustering” by Pepe Berba

“A gentle introduction to HDBSCAN and density-based clustering” by Pepe Berba

“How To Tune HDBSCAN” by Charles Frenzel

Figure 5: The experimental results from ENDLab for (a) single and (c) double mixed layer. Adapted from Ghaemsaidi et al. 2016 [4]. My direct numerical simulation of the (b) single and (d) double mixed layer conditions of the wave tank.

Transmission Coefficient

In order to determine whether internal wave propagation through density staircases is significant for vertical heat transport up from the warm Atlantic layer, it is necessary to quantify how much of an incoming internal wave gets through the staircases. For this, we use the transmission coefficient which is the ratio between the amplitude of the wave that comes out the bottom of the staircase to the amplitude of the original wave that went in the top of the staircase.

Using the Boussinesq equations, Ghaemsaidi et al. 2016 solved for the transmission coefficient and found it depended on two factors: mL, the ratio between the wavelength and thickness of the layer, and θ, the angle at which the wave hits the top of the staircase [4]. As mL gets bigger, the wave gets smaller in comparison to the layer.

Figure 6 shows a heat map of the transmission coefficient Ghaemsaidi et al. 2016 calculated in 4 cases: either 1 or 2 layers and either assuming the fluid is inviscid or viscous [4]. The 1 layer case is exactly what you might intuitively expect; the transmission coefficient drops off as mL gets larger. But, the effect of wave tunnelling can be seen in the 2 layer case where tendrils of high transmission reach across the plots at certain combinations of mL and θ.

Measuring transmission coefficients in numerical experiments

Ghamesaidi et al. 2016 solved for the transmission coefficient analytically [4]. However, doing so becomes much more difficult in more complicated situations, e.g., including both viscosity and the effect of Earth’s rotation. However, adding in such effects is relatively easy in the Dedalus framework. I adapted the code I used above to measure the transmission coefficient for a variety of cases.

To save on computational resources, I collapsed my simulations to 1 spatial dimension. This also allows me to plot an entire simulation in one image as shown in Figure 7. I simulate waves being produced at the top of the domain which propagate downwards through the stratification structure and get absorbed by a sponge layer at the bottom. I run the simulation until it reaches a steady state and then calculate the transmission coefficient using the wavefield enclosed in the “Analyzed domain” box shown in Figure 7.

To calculate the transmission coefficient, I first remove the upward propagating reflected wave from the analyzed domain using a process called complex demodulation, following Mercier et al. 2008 [5]. After that, I take the average amplitude squared above (below) to get a value for the incident (transmitted) wave, then take the ratio to get the transmission coefficient.

Figure 7: An example simulation with 2 layers where mL = 1, θ = 45°, and omega = 0.71. The background stratification profile is shown on the left. The incoming and transmitted waves can be seen separately as indicated by the arrows. The reflected wave can be seen superimposed on the incoming wave. The wavefield enclosed by the “Analyzed domain” box was used to compute a transmission coefficient of 0.28.

Transmission coefficient for multiple layers

Figure 9 shows the transmission coefficient T for many numerical experiments across a certain range of the ratio between wavelength and layer thickness kL at an incident angle of θ = 45°. For the experiments with 1 layer, T decreases monotonically as the layer becomes larger in relation to the wave and matches very well with the analytical prediction. However, for the experiments with multiple layers, there are peaks in transmission at intermediate values of kL. These peaks in transmission come from wave tunnelling and seem to be moving to lower values of kL as the number of layers increases.

Figure 9: The transmission coefficient for many numerical experiments ran in Dedalus across different values of the ratio between wavelength and layer thickness kL for the cases of 1, 2, 3, 4, and 5 layers all with an incident angle of θ = 45°.

Comparing to previous results

Ghaemsaidi et al. 2016 used MATLAB to calculate the analytical transmission coefficients for the cases of 1 and 2 layers [4]. Figure 8 shows a comparison of their results and what I found from numerical experiments for 1 and 2 layers at θ = 45°. The overall shape is quite similar, with the secondary peak in the 2 layer cases corresponding to the main tendril in Figure 6b.

Figure 8: A comparison between the MATLAB code used by Ghaemsaidi et al. 2016 [4] and my simulations using Dedalus for 1 and 2 mixed layers at θ = 45°.

Future work

The work presented here shows that numerical experiments using the Dedalus framework reproduce the results of previous work done using a laboratory wave tank, numerical experiments, or analytical calculations. This implies that such numerical experiments as I’ve shown here could be used to investigate the interactions between internal waves and staircase-like stratification structures in situations that more closely represent the reality of the Arctic Ocean. This would involve using wavelengths and angles of propagation commonly seen in the ocean, which are much larger and shallower than the values I used [4]. This would also mean using stratification structures with values of buoyancy frequency common to the Arctic and interfaces which are much thinner than the layers.

There are many avenues for improving upon the work presented here. For one, I used only weakly-nonlinear, monochromatic waves. However, the internal wave field in the Arctic Ocean contains many different frequencies and nonlinear processes such as harmonic generation can change which frequencies are present and has been shown to modulate internal wave interactions with staircase structures [12]. Bracamontes-Ramirez and Sutherland 2024 showed how transient internal wave packets do not exhibit peaks in transmission as did the steady state in my experiments [1]. To my knowledge, no study has yet explored how the inclusion of both viscosity and Earth’s rotation affects internal wave’s interactions with staircase structures.

References

[1] Bracamontes-Ramirez, J. and B. Sutherland (2024). Transient internal wave excitation of resonant modes in a density staircase. Physical Review Fluids. 9)6):1–21. DOI: 10.1103/PhysRevFluids.9.064801

[2] Eckart, C. (1961). Internal waves in the ocean. Physics of Fluids, 4(7): 791–799. DOI: 10.1063/1.1706408

[3] Fer, I. (2014). Near-Inertial Mixing in the Central Arctic Ocean. Journal of Physical Oceanography, 44 (8): 2031–2049. DOI: 10.1175/JPO-D-13-0133.1

[4] Ghaemsaidi, S. J., H.V. Dosser, , L. Rainville, and T. Peacock (2016). The impact of multiple layering on internal wave transmission. Journal of Fluid Mechanics, 789: 617–629. DOI: 10.1017/jfm.2015.682

[5] Mercier, M., N. Garnier, and T. Dauxois (2008). Reflection and diffraction of internal waves analyzed with the Hilbert transform. Physics of Fluids, 20 (8) 086601. DOI: 10.1063/1.2963136

[6] Shibley, N. and M.L. Timmermans (2019) The Formation of Double-Diffusive Layers in a Weakly Turbulent Environment. Journal of Geophysical Research: Oceans, 124. DOI: 10.1029/2018JC014625

[7] Shibley, N., M.L.Timmermans, J.R. Carpenter, and J.M. Toole (2017). Spatial variability of the Arctic Ocean’s double-diffusive staircase. Journal of Geophysical Research: Oceans, 122 (2):980-994. DOI: 10.1002/2016JC01241

[8] Sutherland, B. R. (2016). Internal wave transmission through a thermohaline staircase. Physical Review Fluids, 1 (1): 013701. DOI: 10.1103/PhysRevFluids.1.013701

[9] Sutherland, B. R. and Yewchuk, K. (2004). Internal wave tunnelling. Journal of Fluid Mechanics, 511: 125–134. DOI: 10.1017/S0022112004009863

[10] Timmermans, M.-L., Toole, J., Krishfield, R., and Winsor, P. (2008). Ice-Tethered Profiler observations of the double-diffusive staircase in the Canada Basin thermocline. Journal of Geophysical Research, 113: 1–10. DOI: 10.1029/2008jc004829

[11] Turner, J. (2010). The melting of ice in the arctic ocean: The influence of Double-Diffusive transport of heat from below. Journal of Physical Oceanography, 40(1):249–256. DOI: 10.1175/2009JPO4279.1

[12] Wunsch, S. (2018). Nonlinear harmonic generation by internal waves in a density staircase. Physical Review Fluids. 3(11):114803. DOI: 10.1103/PhysRevFluids.3.114803