Double‐Diffusive Layer and Meltwater Plume Effects on Ice Face Scalloping in Phase‐Change Simulations

Antarctic ice shelves are losing mass at increasing rates, yet the ice‐ocean interactions that cause significant ice loss are not well understood. A new approach of high‐resolution phase‐change simulations is used to model vertical ice melting into a stratified ocean. The ocean dynamics show complicated interplay between a turbulent buoyant meltwater plume and double‐diffusive layers, while the ice actively melts and changes topography. At room temperatures, the double‐diffusive layer thickness is closely linked to ice scalloping. At lower, more realistic ocean temperatures, the meltwater plume becomes prominent with a laminar‐to‐turbulent transition imprinting an indent on the melting ice. The double‐diffusive layer thickness is consistent with scaling prediction, while the real‐world application demonstrates reasonably good matching of the scaling prediction for some Antarctic regions. Our study is a key first step toward the future use of high‐resolution phase‐change fluid dynamics simulations to better understand Antarctic ice shelves in a changing climate.

Observations reveal rich ocean dynamics (Davis & Nicholls, 2019;Jenkins et al., 2010;Nicholls et al., 2012;Stewart et al., 2019) and ice topography beneath ice shelves (Dutrieux et al., 2014;Schmidt et al., 2023). Some regions show evidence of a buoyant meltwater plume adjacent to the ice (Stanton et al., 2013), which may enhance the melt by entraining warmer and saltier waters toward the ice (Jenkins, 2011). Meltwater discharge from the grounded ice sheet has been linked to ice shelf basal channels, both in Antarctica (Le Brocq et al., 2013) and Greenland (Washam et al., 2019(Washam et al., , 2020. Basal terraces have been observed beneath ice shelves, ranging from smaller than one meter  up to tens of meters (Dutrieux et al., 2014) in the vertical, and up to hundreds of meters in the horizontal. Ice surface scalloping has also been observed beneath some ice shelves (e.g., scallops 0.5 m in width observed by Lawrence et al., 2023), while other ice shelves have been observed to be smooth and flat on their underside (Davis & Nicholls, 2019). Observations beneath George VI (Kimura et al., 2015;Middleton et al., 2022) and Ross (Begeman et al., 2018;Stevens et al., 2020) ice shelves show layering or "steps" in salinity and temperature profiles with depth. Similar steps have been measured near the ocean surface a few hundred meters from an Antarctic glacier tongue, which have been attributed to double-diffusive (DD) layering from sidewall melting (Jacobs et al., 1981). DD processes occur when there are two active scalars with a large difference in molecular diffusivity. Here, the scalars are temperature and salinity, where heat diffuses approximately one hundred times faster than salt. The large diffusivity difference can cause a stable buoyancy profile to evolve into vertical steps, with well-mixed regions broken up by strong gradients in temperature and salinity, resulting in layers (e.g., Radko, 2013). For example, the temperature and salinity steps in the water column beneath George VI Ice Shelf (Kimura et al., 2015). There are two main DD processes that can lead to layering: double-diffusive layers driven by side melting of ice (such as in Jacobs et al., 1981) and double-diffusive convection by ice melting from above (which was attributed to the layering in Kimura et al., 2015;Begeman et al., 2018). However, it is not clear which, or both, of these DD processes occurs beneath different ice shelves, particularly near the grounding line or regions with sloped ice. These complex ocean dynamics (meltwater plume, double-diffusion, ice topography) are expected to have important impacts on grounding line melting (Dutrieux et al., 2013;Marsh et al., 2016;Schmidt et al., 2023). Our focus here is on ice melting from the side to investigate meltwater plume and DD layering effects on ice topography. While vertical ice may not be as common for ice shelf orientation as horizontal or sloped ice, it is still an interesting and relevant case (e.g., the prevalence of vertical and steeply-sloped ice in some regions near Thwaites Eastern and Ross Ice Shelves grounding lines; Schmidt et al., 2023;Lawrence et al., 2023), which has not been considered extensively.
Process-based studies have made progress in examining DD layer, meltwater plume and ice face melting feedbacks. However, these have been primarily limited to laboratory experiments or idealized numerical models in warm ocean temperatures of >0°C, while Antarctic Ocean conditions are generally <0°C. Huppert and Turner (1980) (henceforth HT80) were the first to examine laboratory experiments of a vertical ice face melting into a salinity gradient at room temperature, where the ice scalloped on the same lengthscale as DD layers. More recent laboratory experiments at temperatures of 3-4°C showed that the meltwater plume influenced the melt rate and noted the presence of DD layers, although there was little evidence of the layers contributing to ice scalloping (McConnochie & Kerr, 2016). At similarly cold temperatures (0-6°C), spontaneous thermohaline steps and a meltwater plume were observed in turbulence-resolving numerical simulations of a melting vertical ice face (Gayen et al., 2016). These past simulations approximated the ice-ocean boundary as a flat, static wall, which does not allow for changes in ice shape or turbulence controlled by topographic ice features. Phase-change numerical simulations of a horizontal ice face showed channels in the ice in line with the near-surface imposed current, but this study used only the single scalar of temperature . Regional and large-ocean models are unable to capture the effects of ice face topography and double-diffusive processes (centimeters to tens or even hundreds of meters in vertical scale), and often even the meltwater plume (hundred of meters in the ice-normal direction). Currently, there is inadequate understanding of these fine-scale ice-ocean effects, their influence on the melt rate and parameterization in ocean models. Here, we use a new suite of turbulence-resolving, phase-change numerical simulations to examine the behavior of the vertical ice-ocean interface. For the first time, the simulations include the meltwater plume, double-diffusive layering and changes in ice topography to successfully model key ice face melting effects with applications to ice shelves.

Methodology
Phase-change numerical simulations are used to model a vertical ice face melting into a stratified ocean (Figure 1). The versatile Dedalus numerical solver  is used with a volume-penalization method to model the ice and liquid water phase change (Hester et al., 2020;Hester, Vasil, & Burns, 2021). The phase-change evolution equations call for a very high-resolution grid that makes the simulations computationally expensive. The simulations are therefore restricted to two-dimensions and to laboratory-scale domains of up to one meter. All scales of turbulence are resolved in these direct numerical simulations. The Dedalus numerical solver was previously used for a melting iceberg application, where it was successfully validated by replicating laboratory experiments (Hester, McConnochie, et al., 2021). A very similar set-up of the Dedalus numerical solver is used here and so the reader is referred to Hester, McConnochie, et al. (2021) for further information about the governing equations, solver set-up, phase-change equations and other numerical details. The only notable difference here is that a trigonomic discretion method is used in both the horizontal and vertical directions to have closed boundaries on all four sides of the domain.
The simulated rectangular domain has length, L, and height, H, with gravity, g, pointing downwards in the z-direction (Table 1). The four boundaries are impermeable and insulating in temperature and salinity. The majority of the domain is initialized as liquid water, with a 2-4 cm wide ice block on the left-hand side spanning the entire vertical domain. The ocean temperature is initialized as homogeneous at the chosen reference temperature,  Table 1). Snapshots are taken 15 min into the laboratory experiment and 5 min into the simulation. The ice on the left-hand side (black in (a) and gray in (b, c)) shows scalloping with the double-diffusive layers.
Case Note. The melt rate, m, has been time-averaged after the first minute of simulation time and the uncertainty range is the standard deviation. For Cases L1-L3, h is the layer thickness measured in the turbulent plume region.

Table 1
Key Parameters for the Simulations and Previous Observational Studies, Where J81 Is Jacobs et al. (1981), K15 Is Kimura et al. (2015) and B18 Is Begeman et al. (2018) T ∞ . The ocean salinity is initialized with a background linear salinity gradient, which is described in terms of a reference salinity, S ∞ (the salinity at mid-height), and the buoyancy frequency, = √ (− ∕ 0) ∕ , where ρ is the density and ρ 0 is the reference density. A nonlinear equation of state relates density to the salinity and temperature (McDougall et al., 2009). The initial temperature of the ice is chosen to be near the expected freezing point for the ocean reference temperature and salinity. The thermal driving is ΔT = T ∞ − T w , where T w is the approximated temperature at the ice face calculated on the freezing point depression curve with S ∞ .
The molecular viscosity is the realistic value of ν = 1.3 × 10 −6 m 2 s −1 . The molecular diffusivities are κ T = 1.3 × 10 −6 m 2 s −1 for heat and κ S = 2.6 × 10 −8 m 2 s −1 for salt, resulting in a Lewis number of Le = κ T /κ S = 50. The real ocean has κ T ≈ 1.4 × 10 −7 m 2 s −1 and κ S ≈ 7 × 10 −10 m 2 s −1 , giving a Lewis number of Le = 180. The salt diffusivity is larger in the model than in the real ocean, as the simulations would require extremely high grid resolution for realistic salt diffusivities. However, if only the model salt diffusivity is increased, then there may not be a large enough difference between salt and heat diffusivities to allow for realistic DD processes. To overcome this, the model thermal diffusivity is increased by an order of magnitude to attain Le = 50. The increased diffusivities may lead to higher heat and salt fluxes, and an enhanced melt rate. For example, Case H3 appears to melt up to three times faster than the laboratory experiment (comparing Figure 1 panels a and b). Similar elevated salt diffusivities are used in past fine-scale numerical studies of ice-ocean interactions (Gayen et al., 2016;Hester, McConnochie, et al., 2021;Mondal et al., 2019).
Scaling arguments from HT80 predict the DD layer thickness, η, to be A scaling constant, c = h/η, is introduced to relate the measured height, h, to the scaling prediction in Equation 1. For the room temperature laboratory experiments, HT80 found c = 0.65. The onset of motion in the layers (also referred to as cellular convection) occurs when the Rayleigh number, is larger than O(10 4 ) − O(10 5 ) (Chen et al., 1971).
Four simulations (H1-H4) are conducted at high temperatures to model HT80 experiments and validate the simulations. Three simulations (L1-L3) are conducted at low temperatures to more closely model the real-world scenario. Most simulations have 2,048 × 2,048 grid cells. The exceptions are Case L2 (2,048 × 4,096) and Case L3 (4,096 × 8,192) which require larger domains for the meltwater plume to become turbulent. As the domain height is restricted by computational resources to around a meter, the salinity stratification is chosen so that at least three to four layers are expected in the domain. This means that the salinity stratification is stronger than is likely for sub-ice regions, which is thought to only impact setting the layer lengthscale and is an approach consistent with HT80. The exception is Case L3, which is similar in magnitude to some stratifications measured beneath Thwaites Eastern Ice Shelf . Simulations are run for between 2.8 and 23 min, in order to see the growth of the meltwater plume and the double-diffusive layers.

Results and Discussion
The ice begins as a rectangular block that melts into the stratified ocean, cooling and freshening adjacent waters.
Next to the ice, the freshening effect dominates the density change over the cooling effect, leading to buoyant boundary layer flow indicative of a meltwater plume. In the lower temperature cases, the plume rises upwards against the ice, where it begins as laminar before becoming turbulent (at the height of 30 cm in Figure 2 and Movie S2). In the higher temperature cases, the meltwater plume is present but less obvious as the dynamics are quickly dominated by strong double-diffusive layers (Figure 1 and Movie S1). The double-diffusive layers are clearest in the higher temperature cases (Figure 1) but are also present in lower temperature cases as interleaving layers moving horizontally outwards from the ice (red box in Figure 2). The mechanism causing double-diffusive layers will now be outlined. After the meltwater plume initializes, the ice continues to melt and the relatively warm and saline water in the ambient far-field starts transferring heat and salt toward the meltwater plume. The large difference in diffusivities means that the heat from the ambient water is diffused toward the ice at a faster rate than salt. The water parcels neighboring the meltwater plume are cooled faster than they are freshened, making them relatively dense. These water parcels begin to sink, but the descent is limited by the stabilizing density (salinity) gradient. Once the water parcels reach a depth of equal density class, they are redirected to move horizontally outwards from the ice, which forms the distinct double-diffusive layers (Figures 1 and 2; Movies S1 and S2). The double-diffusive layers also appear as steps in the salinity profile, which are well-mixed regions followed by sharp gradients (Case L3 in Figure 2d shows sharp salinity gradients of dS/dz ≈ 0.03 PSU cm −1 ). In the temperature profile, the DD layers appear as sharp changes in gradient (Case L3 in Figure 2d shows sharp temperature gradients of dT/dz ≈ 0.3°C cm −1 ). Over time, the DD layers continue to expand further outwards as the layer circulation draws in ambient water further from the ice.
After the DD layers form, a secondary instability of salt finger plumes begins within the layers (plumes clearly shown in Figure 1b). Salt fingering is another type of DD instability that forms when the overall density gradient is stable, but there is warmer and saltier water of lower density overlying colder, fresher water of higher density (Radko, 2013). In our simulations, the incoming flow within a layer is relatively warm and salty compared to the underlying outgoing flow, while the locally stable density gradient is maintained. Salt fingering keeps the DD layers well-mixed between the sharp temperature and salinity gradients (Figures 1 and 2).
At higher temperatures, the ice undergoes significant scalloping with vertical (z-direction) topographic lengthscale matching the DD layer thickness, consistent with the HT80 experiments (approximately 4 cm in Figure 1). The melt rate is enhanced at the trough of the scallops (Figure 1d), which continues deepening the scallops over time. This behavior is similar to laboratory experiments of near-horizontal ice melting into a free-stream flow, where scallops also spontaneously formed with melt rates highest in the trough (Bushuk et al., 2019). In our simulations at lower temperatures, there is no clear signal of DD layers directly impacting the ice topography (Figure 2), as also found in experiments by McConnochie and Kerr (2016). The lower temperature cases instead have an indent in the ice at the height where the meltwater plume transitions from laminar to turbulent. For example, in Case L3, the plume becomes turbulent at z ≈ 30 cm and the subtle indent (5 mm deep) follows at z ≈ 35 cm (Figure 2a). Compared to the laminar region, the turbulent plume has an increased upward velocity that draws in the far-field water at the transition point. The far-field water is relatively warmer and saltier than the cooled and freshened water close to the ice-ocean boundary, resulting in a locally increased melt rate (Figure 2b). At higher temperatures, the DD layers do not evolve in thickness over time, allowing them to imprint on the ice face (Figure 3a and Figure S1 in Supporting Information S1). At lower temperatures, the layers evolve in thickness, with new layers forming and older layers merging or disappearing (Figure 3b and Figure S1 in Supporting Information S1). The layer modulation may be due to the increased relative strength of the turbulent meltwater plume at lower temperatures, compared to the strength of the DD layers (see Movies S1 and S2 for comparison of plume and DD layer strength). This comparatively strong plume may broaden, oscillate and move around the circulation pathways near the ice.
The DD layer thickness is broadly consistent with the scaling (Equation 1; Figure 3c), meaning that the layer thickness can be reasonably well predicted using only the temperature, salinity and stratification of the far-field ocean. However, the lower temperature cases exhibited some difference in layer behavior between the laminar and turbulent plume regions, with the latter having thicker layers that progress outwards from the ice more quickly (Figure 2a). For this reason, the layer thickness in the laminar and turbulent plume regions have been included separately on Figure 3c. We speculate that the thicker layers are due to the turbulent meltwater plume entraining with neighboring waters, which enhances the heat and salt mixing above molecular diffusion values in this region. The added mixing can be thought of as a turbulent diffusivity that equally affects the temperature and salinity, thereby resulting in a lower effective Lewis number, weaker double-diffusive processes, and thicker layers.

Ocean Application
The real-world application is Antarctic ice shelves melting into a stratified ocean. There are three main types of instabilities that have been used to explain observed layer thickness in these regions: (a) double-diffusive layers from the side melting of ice which is the focus of the present study, (b) double-diffusive convection due to ice melting from above (Radko, 2013), and (c) stable stratification plus mixing from an external source such as internal wave breaking (with no DD processes involved; Park et al., 1994). It is difficult to untangle these three mechanisms in different Antarctic regions, hence we rely on a comparison with the scaling prediction to indicate whether DD layering (mechanism (a)) is a candidate for the observed layers.
The DD layer (mechanism (a)) scaling does well predicting the layer thickness for observations taken by Jacobs et al. (1981) near the ocean surface adjacent to an Antarctic glacier tongue, as also noted by the authors of that study (Table 1). The profiles (reproduced in Figure 3d) show stable salinity and temperature gradients throughout most of the water column. As double-diffusive convection (mechanism (b)) relies on the formation of an unstable temperature gradient, it is unlikely to explain the Jacobs et al. (1981) case. A more recent observational study in this region did not find additional evidence of layering, but did outline that mechanism (c) may also explain the previously observed layers, noting strong tidal motions moving across bottom topography as a possible source of water column mixing (Stevens et al., 2014). Observations of layers beneath George VI (Kimura et al., 2015) and Ross (Begeman et al., 2018) ice shelves both showed a relatively quiescent ocean with low mean currents and tidal motions, making mechanism (c) unlikely for this data. The George VI Ice Shelf observations do not match that well with the DD layer (mechanism (a)) scaling (Table 1). Instead, and as noted by Kimura et al. (2015), double-diffusive convection (mechanism (b)) is a reasonable explanation for the layer observations beneath this relatively horizontally flat region of the ice shelf. Double-diffusive convection creates a well-mixed layer immediately beneath horizontally flat ice (Middleton et al., 2021) and can contribute to further layers being created beneath (Rosevear et al., 2021), which could result in the steps in the temperature and salinity profiles beneath George VI Ice Shelf. The steps beneath George VI Ice Shelf were observed within 10-20 m of the ice rather that further down in the water column, which is also consistent with double-diffusive convection layers forming from above (mechanism (b)). For Ross Ice Shelf, the observations by Begeman et al. (2018) were taken near the grounding line, a region where the ice shelf may be more sloped. The DD layer (mechanism (a)) scaling is somewhat reasonable in comparison to the observed layer thickness (Table 1). But overall, the result remains rather inconclusive, as Begeman et al. (2018) find that the water column is likely unstable to double-diffusive convection (mechanism (b)). It is worth noting that double-diffusive convection by ice melting from above (mechanism (b)) is a relatively slow process, while DD layering by side melting (mechanism (a)) has layers at depth evolving outwards at the same rate. The case of a vertical ice face melting into a stratified ocean (mechanism (a)) may then be of importance beneath other Antarctic ice shelves, in particular when considering that the basal topography of ice shelves can have channels, crevasses, or terraces from meters to hundreds of meters in vertical scale (Dutrieux et al., 2014;Lawrence et al., 2023;Schmidt et al., 2023).
Our simulations show scalloping of the ice at higher temperatures, but not at lower temperatures relevant to the ocean. However, double-diffusive layers in the ocean may still be able to impact ice topography. The ocean is at significantly larger length-and time-scales than modeled here and subtle changes in the topography may be enhanced over larger distances and longer times. On the other hand, the layers may modulate and move around in time, as noted in the lower temperature simulations, which may mean that an imprint is not observed on the ice even at large length-and time-scales. The George VI Ice Shelf observations showed that the double-diffusive layers changed significantly over only an hour or so (Kimura et al., 2015). Even if the double-diffusive layers do not imprint on the ice topography, they still have a significant effect on the stratification and the transport of heat and salt from the far-field toward the ice, which will likely impact the melt rate. Finally, we note that our primary focus has been on Antarctic ice shelves, but the mechanisms described in this study could also be applicable to ice shelves in the Arctic.

Conclusions
For the first time, high-resolution phase-change numerical simulations are used to model ice-ocean interactions with a vertical ice orientation, while resolving double-diffusive processes and changing ice topography. Vertical 10.1029/2023GL104396 8 of 9 ice melting produces double-diffusive layers at both room temperatures and lower temperatures more representative of Antarctica. A buoyant meltwater plume is particularly clear at lower temperature, with a laminar to turbulent transition resulting in an indent in the ice. Ice scalloping on the topographic scale of the double-diffusive layers occurred at room temperatures, but not in lower temperature cases. However, we speculate that the real-world ice-ocean application may yet see ice scalloping from double-diffusive layers, particularly across larger distances and longer times than were accessible in our numerical simulations.
Our study demonstrates the emerging possibilities of high-resolution phase-change numerical simulations to model key ice-ocean interactions. Future increases in computing power will allow more realistic length-and time-scales. Further simulations could include additional processes such as across ice-face currents (constant or tidal) and glacial outflow, while sloped ice would be relevant to many ice shelf regions. In addition to a numerical approach, laboratory experiments remain extremely useful in understanding the dynamics and observational studies continue making important, new measurements of the real system. Improved knowledge of ice-ocean interactions from all these approaches will update parameterizations in regional and large scale ocean models, to better understand ice shelf and ocean effects in a changing climate.

Data Availability Statement
The Dedalus solver (version 2) was used for all the numerical simulations. The Dedalus solver is developed openly and available at http://dedalus-project.org. Model output is preserved at https://doi.org/10.5281/zenodo.8175219.