Large-eddy simulations of the ice-shelf–ocean boundary layer near the ice front of Nansen Ice Shelf, Antarctica

Ice melting beneath Antarctic ice shelves is caused by heat transfer through the ice-shelf–ocean boundary layer (IOBL). However, our understanding of the fluid dynamics and thermohaline physics of the IOBL flow is poor. In this study, we utilize a large-eddy simulation (LES) model to investigate ocean dynamics and the role of turbulence within the IOBL flow near the ice front. To simulate the varying turbulence intensities, we imposed different theoretical profiles of the velocity. Far-field ocean conditions for the melting at the ice-shelf base and freezing at the sea surface were derived based on in situ observations of temperature and salinity near the ice front of the Nansen Ice Shelf. In terms of overturning features near the ice front, we validated the LES simulation results by comparing them with the in situ observational data. In the comparison of the velocity profiles to shipborne lowered acoustic Doppler current profiler (LADCP) data, the LES-derived strength of the overturning cells is similar to that obtained from the observational data. Moreover, the vertical distribution of the simulated temperature and salinity, which were mainly determined by the positively buoyant meltwater and sea-ice formation, was also comparable to that of the observations. We conclude that the IOBL flow near the ice front and its contribution to the ocean dynamics can be realistically resolved using our proposed method. Based on validated 3D-LES results, we revealed that the main forces of ocean dynamics near the ice front are driven by positively buoyant meltwater, concentrated salinity at the sea surface, and outflowing momentum of the sub-ice-shelf plume. Moreover, in the strong-turbulence case, distinct features such as a higher basal melt rate (0.153 m yr⁻¹), weak upwelling of the positively buoyant ice-shelf water, and a higher sea-ice formation were observed, suggesting a relatively high speed current within the IOBL because of highly turbulent mixing. The findings of this study will contribute toward a deeper understanding of the complex IOBL-flow physics and its impact on the ocean dynamics near the ice front.

1 Introduction

The Antarctic Ice Sheet (AIS) is buttressed by floating extensions of land ice, called ice shelves (Rignot et al., 2013), which primarily control the AIS mass balance by hindering the flow of inland ice into the ocean (Holland et al., 2020). Iceberg calving and basal melting are the key factors that can weaken the buttressing capacity of ice shelves (Liu et al., 2015; Smith et al., 2020).

The sub-ice-shelf oceanic environment can be divided into two classes, namely the “cold-water cavity” (e.g., the Larsen C, Ross, Filchner–Ronne, and Amery ice shelves) and the “warm-water cavity” (e.g., the Getz and Totten ice shelves), depending on the amount of basal melting as well as ocean conditions (Gwyther et al., 2016; Joughin et al., 2012). In the cold-water cavity, shear force generated by the tides and sinking force of dense waters generated by brine rejection (e.g., high-salinity shelf water, HSSW) are the driving forces that cause basal melting (Davis and Nicholls, 2019; Yoon et al., 2020). In contrast, the intrusion of circumpolar deep water and melt-driven circulation near the grounding line are major causes of strong basal melting in the warm-water cavities (Holland et al., 2020; Jacobs et al., 1992). Investigating the ice-shelf–ocean boundary layer (IOBL), which is the boundary layer (meters to tens of meters) right beneath the ice shelf, is a complex problem because turbulent mixing and the stably stratified ocean layer generated by basal melting both influence the IOBL characteristics (Begeman et al., 2018; Garabato et al., 2017). Various observational studies beneath ice shelves of Antarctica have been performed to observe the thermohaline characteristics and structures of the IOBLs as well as identify the ocean conditions beneath the ice shelves (Jenkins et al., 2010; Kimura et al., 2015). In the Larsen C ice shelf, which is a cold-water cavity, a well-mixed boundary layer (20–30 m) was observed in both temperature and salinity, induced by a strong tidal forcing and a weak stratification. A moderate melt rate (1.9 m yr⁻¹) was observed, despite the low thermal driving due to the observed shear-driven turbulence (Davis and Nicholls, 2019). Similar melt rates for weak stratification were also observed beneath the Fimbul and Ross ice shelves (Arzeno et al., 2014; Hattermann et al., 2012). Moreover, the heat entrainment is prevented near the grounding line of the Ross Ice Shelf because of the strong density gradient between the meltwater and the denser continental shelf water in this region, yielding a low melt rate (Begeman et al., 2018).

Positively buoyant sub-ice-shelf plumes that are created near the grounding line and rise up along the ice-shelf bottom can affect the stratification and heat entrainment within the IOBL because the upper part of sub-ice-shelf plumes is the IOBL (Hewitt, 2020; Holland and Jenkins, 1999). Since the sub-ice-shelf plume modulates the IOBL characteristics with a strong turbulent mixing, the basal melt rate can be increased by weakening the stratification and enhancing the entrainment of the seawater from the outer region. As the sub-ice-shelf plume generated by the HSSW in the cold-water cavities gets closer to the ice front, it flows outward, losing its buoyancy (Smethie and Jacobs, 2005).

Observational efforts of meltwater behavior and ocean circulation near the frontal region of the ice shelf demonstrate various mechanisms at different locations around Antarctica. In the frontal region of the Pine Island Ice Shelf, Garabato et al. (2017) revealed that the ascent of the meltwater outflow causes vigorous lateral export, affecting the settling of meltwater at depth. The intrusion of relatively warm surface waters and high basal melting near the ice-shelf front was observed in the Ross Ice Shelf (Horgan et al., 2011). Moreover, Malyarenko et al. (2019) suggested the existence of a “wedge” of fresher water in the western Ross Sea, and it is formed from meltwater near the ice-shelf front. However, understanding of meltwater generation and behavior near the frontal region of the ice shelf is still poor, because the observation of sub-ice-shelf environments is quite limited.

In contrast to the Reynolds-averaged Navier–Stokes model which provides a time-averaged feature, high-resolution turbulence models such as the large-eddy simulation (LES) or direct numerical simulation (Gayen et al., 2016; McConnochie and Kerr, 2018; Mondal et al., 2019; Vreugdenhil and Taylor, 2019) can resolve the detailed eddy structures and quantify the momentum and heat fluxes near the wall boundary. Therefore, the LES for the IOBL flow is required to analyze the detailed IOBL structure and reveal fundamental information about the basal melting phenomenon occurring below the ice shelves (Jenkins, 2016; Vreugdenhil and Taylor, 2019). The LES has been successfully applied to investigate the role of turbulence in the ice–ocean interactions in Greenland (Carsey and Garwood, 1993; Denbo and Skyllingstad, 1996) and the Arctic Sea (Glendening, 1995; Skyllingstad and Denbo, 2001; Matsumura and Ohshima, 2015; Ramudu et al., 2018; Li et al., 2018). However, studies on the application of LES to the IOBL under a sub-ice-shelf environment are limited due to lack of observations (Dinniman et al., 2016).

In this study, we performed LES experiments for the IOBL and ocean flow with neutrally buoyant sub-ice-shelf plumes near the ice front. To include the thermohaline effect by sea-ice formation at the sea surface and basal melting at the ice-shelf base, surface fluxes in both temperature and salinity were used. The boundary conditions used in LES experiments were derived based on the in situ observation data – namely the conductivity–temperature–depth (CTD), lowered acoustic Doppler current profiler (LADCP) data, and automatic weather station (AWS) data – collected in front of the Nansen Ice Shelf (NIS), which is a cold-water cavity, during the Antarctic expedition led by the Korea Polar Research Institute in January and February 2017. One of the main objectives of this study was to investigate the neutrally buoyant sub-ice-shelf plume's impact on basal melting. Therefore, the target parameter in the LES experiments was turbulence intensity within the IOBL because the sub-ice-shelf plume and its heat entrainment are related to basal melting via turbulent shear. Another objective was to validate our proposed methodology (domain configuration and boundary conditions) and the oceanographic properties simulated using the LES. Using the validated three-dimensional LES outputs, we quantified the distribution of melting at the ice-shelf bottom as well as the associated factors such as turbulent characteristics and flux changes within the IOBL.

Section 2 of this paper presents the governing equations (i.e., the Navier–Stokes equation and liquidus condition) for the oceanic flow with melting and freezing effects as well as a detailed explanation of the simulations. Section 3 presents the analysis of the LES simulation results to determine the IOBL characteristics – flow velocity, potential temperature, salinity, fluxes, and turbulence statistics. The major findings, future works, and implications of this study are summarized in Sects. 4 and 5.

2 Methodology

2.1 CTD and LADCP observations

During the hydrographic survey by the ice-breaking research vessel ARAON operated by the Korea Polar Research Institute, full-depth CTD and LADCP profiles were obtained in 1 h intervals from 14 February 2017 at 13:00 UTC to 15 February 2017 at 11:55 UTC at a single site in the Terra Nova Bay in front of the NIS. This survey was conducted to examine the vertical structures of the sub-ice-shelf plume with temporal variations (grey lines in Fig. 5). The exact location of the observations was 75.008^∘ S, 163.617^∘ E (Fig. 1b), located approximately 1 km away from the ice front. Because the LADCP data were incorrectly observed at the final casting, 24 CTD profiles and 23 LADCP profiles were used in this study. The CTD data were processed following the SBE recommended procedure (Sea-Bird Electronics, Inc., Bellevue, Washington, USA, 2014), and the LADCP data were processed using the methods introduced in Thurnherr (2014). This location for the observations is the coastal polynya region where strong katabatic winds cause extreme heat loss in the ocean, leading to sea-ice formation. Atmospheric properties such as wind speed, temperature, sensible heat, and sea-ice formation were obtained using the AWS instrument on ARAON and are listed in Table 1. Based on the wind speeds and air temperatures in the AWS data acquired on ARAON, we calculated the sensible heat as well as the amount of sea-ice formation (Thompson et al., 2020). The detailed shipboard information and processing methods for the hydrographic data are described in Yoon et al. (2020).

Table 1 List of model parameters and constants.

^a Heorton et al. (2017). ^b Based on friction velocity (0.168 m s⁻¹) and thermal driving (refer to Vreugdenhil and Taylor, 2019). ^c Smooth ice with melting case (C_d=0.001) in Gwyther et al. (2016). ^d Stevens et al. (2017). ^e Thompson et al. (2020).

Figure 1 Region of interest and simulation domain configuration. (a) Map of Antarctica. The red box shows the study area, Terra Nova Bay, in the western Ross Sea. (b) Region of interest where the CTD and LADCP surveys were conducted in the 2016–2017 shipboard survey (satellite image was obtained from © Google Earth, 2020). (c) Simulation domain and boundary conditions.

2.2 Numerical model

To simulate the oceanic flow incorporating the effects of sea-ice formation and basal melting, the parallelized large-eddy simulation model (PALM, version 6–r4536) developed by Leibniz University was employed (Noh et al., 2009; Raasch and Schröter, 2001). This model solves the non-hydrostatic, Boussinesq-approximated, filtered Navier–Stokes equations with buoyancy force, Coriolis force, and subgrid-scale (SGS) turbulent closure. The Boussinesq approximation can be applied to flows with negligible density variation. Furthermore, in the time integration, the time-difference formulas were computed using the third-order Runge–Kutta method. The fifth-order upwind scheme was used to solve the flow advection (Wicker and Skamarock, 2002). The pressure was modeled using a Poisson equation, while the mass, momentum, potential temperature, and salinity conservations were governed by Eqs. (1)–(4) (Einstein summation convention), respectively.

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_k}}{\partial x_k}}=\mathrm{0},{\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_{\mathrm{i}}}}{\partial t}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_{\mathrm{i}}}\stackrel{\mathrm{‾}}{u_k}}{\partial x_k}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\pi }}^{\ast }}}{\partial x_k}}-{\mathit{\epsilon }}_{ijk}{f}_j\stackrel{\mathrm{‾}}{u_k}+{\mathit{\epsilon }}_{i\mathrm{3}k}{f}_{\mathrm{3}}\stackrel{\mathrm{‾}}{u_{g,k}}\\ {\displaystyle }& {\displaystyle }+g{\displaystyle \frac{\left({\mathit{\rho }}_{\mathit{\theta }}-<{\mathit{\rho }}_{\mathit{\theta }}>\right)}{<{\mathit{\rho }}_{\mathit{\theta }}>}}{\mathit{\delta }}_{i\mathrm{3}}-{\displaystyle \frac{\partial {\mathit{\tau }}_{ki}^{\mathrm{r}}}{\partial x_k}},\\ \text{(2)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{ki}^{\mathrm{r}}={\mathit{\tau }}_{ki}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\mathit{\tau }}_{jj}{\mathit{\delta }}_{ki},\stackrel{\mathrm{‾}}{{\mathit{\pi }}^{\ast }}=\stackrel{\mathrm{‾}}{p^{\ast }}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\mathit{\tau }}_{jj}{\mathit{\delta }}_{ki},\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\theta }}}{\partial t}}=-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_k}\stackrel{\mathrm{‾}}{\mathit{\theta }}}{\partial x_k}}+{\displaystyle \frac{\partial H_k}{\partial x_k}}+Q_{\mathit{\theta }},\text{(4)}& {\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{S}}{\partial t}}=-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_k}\stackrel{\mathrm{‾}}{S}}{\partial x_k}}+{\displaystyle \frac{\partial S_k}{\partial x_k}}+Q_{\mathrm{s}},\end{array}$$

where u_k is the flow velocity, ρ is the seawater density, π^∗ is the dynamic pressure, ε_ijk is the Levi–Civita symbol, f is the Coriolis force for 75S (−1.41 × 10⁻⁴ s⁻¹), δ_ij is the Kronecker delta function, g is the gravitational acceleration, ρ_θ is the potential density, T is the absolute temperature, ν is the dynamic viscosity, ${\mathit{\tau }}_{kj}^{\mathrm{r}}$ is the Reynolds stress, θ is the potential temperature, p is the hydrostatic water pressure, p₀ is the reference pressure, and S is the salinity (Jackett et al., 2006). Additionally, Q_θ and Q_s are the external forcing of the source and sink terms – T and S, respectively. The overbars indicate that the values have been filtered over the grid volume. Combining these equations, the SGS turbulent kinetic energy (e) equation can be derived as follows:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{e}}{\partial t}}& =-\stackrel{\mathrm{‾}}{u_k}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{e}}{\partial x_k}}-{\mathit{\tau }}_{ki}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_{\mathrm{i}}}}{\partial x_k}}\\ & +{\displaystyle \frac{g}{{\mathit{\theta }}_{\mathrm{0}}}}\mathit{\{}\stackrel{\mathrm{‾}}{u_{\mathrm{3}}^{\prime }{\mathit{\rho }}_{\mathit{\theta }}^{\prime }}-{\displaystyle \frac{\partial }{\partial x_k}}\stackrel{\mathrm{‾}}{u_k^{\prime }\left(e+{\displaystyle \frac{p^{\prime }}{{\mathit{\rho }}_{\mathrm{0}}}}\right)}\mathit{\}}-\mathit{\epsilon },\end{array}\end{array}$$

where $=\frac{{\stackrel{\mathrm{‾}}{u^{\prime }}}_{\mathrm{i}}{\stackrel{\mathrm{‾}}{u^{\prime }}}_{\mathrm{i}}}{\mathrm{2}}$ and ε is the SGS dissipation rate.

The SGS stresses (τ_ki, H_k, and S_k) for the momentum, potential temperature, and salinity are parameterized as follows (Maronga et al., 2015):

$$\begin{array}{ll}\text{(6)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{ki}=\stackrel{\mathrm{‾}}{u_k{u}_{\mathrm{i}}}-\stackrel{\mathrm{‾}}{u_k}\stackrel{\mathrm{‾}}{u_{\mathrm{i}}}=-K_{\mathrm{m}}\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_{\mathrm{i}}}}{\partial x_k}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_k}}{\partial x_{\mathrm{i}}}}\right)+{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\mathit{\delta }}_{ik}\stackrel{\mathrm{‾}}{e},{\displaystyle }& {\displaystyle }{H}_k=\stackrel{\mathrm{‾}}{u_k\mathit{\theta }}-\stackrel{\mathrm{‾}}{u_k}\stackrel{\mathrm{‾}}{\mathit{\theta }}=\stackrel{\mathrm{‾}}{u_k^{\prime }\stackrel{\mathrm{´}}{\mathit{\theta }}}=-K_{\mathrm{h}}\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\theta }}}{\partial x_k}}\right),\\ \text{(7)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{h}}=\left(\mathrm{1}+\mathrm{2}{\displaystyle \frac{l}{\mathrm{\Delta }}}\right)K_{\mathrm{m}},\text{(8)}& {\displaystyle }& {\displaystyle }{S}_k=\stackrel{\mathrm{‾}}{u_{k}S}-\stackrel{\mathrm{‾}}{u_k}\stackrel{\mathrm{‾}}{S}=\stackrel{\mathrm{‾}}{u_k^{\prime }{S}^{\prime }}=-K_{\mathrm{h}}\left({\displaystyle \frac{\partial S}{\partial x_k}}\right),\text{(9)}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial x_k}}\left[\stackrel{\mathrm{‾}}{u_k^{\prime }(e+{\displaystyle \frac{p^{\prime }}{{\mathit{\rho }}_{\mathrm{0}}}})}\right]=-{\displaystyle \frac{\partial }{\partial x_k}}{\mathit{\upsilon }}_e{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{e}}{\partial x_k}},{\mathit{\upsilon }}_e=\mathrm{2}{K}_{\mathrm{m}},{\displaystyle }& {\displaystyle }\mathit{\epsilon }=C_{\mathit{\epsilon }}{\displaystyle \frac{e^{\mathrm{3}/\mathrm{2}}}{l}},C_{\mathit{\epsilon }}=\mathrm{0.19}+\mathrm{0.74}l,\\ {\displaystyle }& {\displaystyle }{K}_{\mathrm{m}}=C_{\mathrm{m}}l\sqrt{\stackrel{\mathrm{‾}}{e}}\text{ with empirical value }{C}_{\mathrm{m}}=\text{constant}=\mathrm{0.1},\\ {\displaystyle }& {\displaystyle }\text{ and }l=\min \left(\mathrm{1.8}z,\mathrm{\Delta },\mathrm{0.76}\sqrt{\stackrel{\mathrm{‾}}{e}}{\left[{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathit{\theta },\mathrm{0}}}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathit{\theta }}}}{\partial z}}\right]}^{-\frac{\mathrm{1}}{\mathrm{2}}}\right),\end{array}$$

where K_m and K_h are the eddy diffusivities for momentum and heat; l is the turbulent mixing length which depends on height z (distance from the wall), grid spacing, and stratification; Δ= ($\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z{)}^{(\mathrm{1}/\mathrm{3})}$ is the length scale of the filter; and ρ_θ is the potential density. (the variables with prime are SGS variables). Because this SGS model and the coefficients are designed for a stably stratified boundary layer flow, it is suitable for resolving a stratified ocean flow with melting or freezing effect. However, the SGS fluxes can be incorrectly determined at a region where the flow becomes laminar, as this model assumes only a turbulent flow.

It is necessary to determine the ambient variables (θ_f and S_f) away from the ice-shelf base and the interfacial variables (θ_b and S_b) near the ice-shelf–ocean boundary to resolve the thermal and salinity changes caused by the freezing effect at the sea surface or the basal melting effect at the ice-shelf–ocean boundary. Herein, we determined the ambient variables based on in situ CTD observations obtained 1 km away from the ice front (Fig. 1b). To obtain the interfacial variables, we solved the conservation equations of heat and salt, along with the liquidus condition and turbulent transfer coefficients for heat and salt (Beckmann and Goosse, 2003; Vreugdenhil and Taylor, 2019).

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{c}_{\mathrm{w}}{\mathit{\rho }}_{\mathrm{w}}{\mathrm{\Gamma }}_{\mathit{\theta }}{u}_{\ast }({\mathit{\theta }}_{\mathrm{f}}-{\mathit{\theta }}_{\mathrm{b}})={\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{i}}m,\text{(11)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\mathrm{w}}{\mathrm{\Gamma }}_{\mathrm{S}}{u}_{\ast }(S_{\mathrm{f}}-S_{\mathrm{b}})={\mathit{\rho }}_{\mathrm{i}}{S}_{\mathrm{b}}m,\text{(12)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{b}}={\mathit{\lambda }}_{\mathrm{1}}{S}_{\mathrm{b}}+{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}P,\text{(13)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{b}}={\displaystyle \frac{\begin{array}{c}-\left({\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}P-{\mathit{\theta }}_{\mathrm{f}}+\frac{L_{\mathrm{i}}{\mathrm{\Gamma }}_{\mathrm{s}}}{c_{\mathrm{w}{\mathrm{\Gamma }}_{\mathit{\theta }}}}\right)\\ +\sqrt{{\left({\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}P-{\mathit{\theta }}_{\mathrm{f}}+\frac{L_{\mathrm{i}}{\mathrm{\Gamma }}_{\mathrm{s}}}{c_{\mathrm{w}{\mathrm{\Gamma }}_{\mathit{\theta }}}}\right)}^{\mathrm{2}}-\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}\left(\frac{L_{\mathrm{i}}{\mathrm{\Gamma }}_{\mathrm{s}}}{c_{\mathrm{w}{\mathrm{\Gamma }}_{\mathit{\theta }}}}\right)S_{\mathrm{f}}}\end{array}}{-\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}},\end{array}$$

where m is the melt rate at the ice-shelf base or the freezing rate at the sea surface; the subscripts w and i refer to the parameters for water and ice, respectively. The parameters values are listed in Table 1. The friction velocity at the ice-shelf base was calculated from the simulated velocity field. We used 0.026 m s⁻¹ as the friction velocity for the calculation of thermal and salinity fluxes induced by the sea-ice formation, although the effect of wind stress on the momentum was excluded to focus on the relationship between the sub-ice-shelf plume and the development of overturning cell (counterclockwise direction cell observed in Fig. 5a).

The fluxes for temperature and salinity, q_θ∗ and q_S∗, at the ice-shelf bottom were formulated by Monin–Obukhov similarity and interfacial values, θ_b and S_b, obtained by the resulting equation, Eq. (13) (McPhee et al., 1987; Ramudu et al., 2018):

$$\begin{array}{}\text{(14)}& {\displaystyle }{q}_{\mathit{\theta }\ast }={\mathrm{\Gamma }}_{\mathit{\theta }}[\mathit{\theta }\left(z_{\mathrm{1}}\right)-{\mathit{\theta }}_{\mathrm{b}}]u_{\ast },\text{(15)}& {\displaystyle }{q}_{\mathrm{S}\ast }={\mathrm{\Gamma }}_{\mathrm{S}}[S\left(z_{\mathrm{1}}\right)-S_{\mathrm{b}}]u_{\ast },\end{array}$$

where z₁ is first node from the boundary; u_∗ is the friction velocity, which is calculated using the logarithmic law of the wall with the velocity at the first node and surface roughness (z₀); and Γ_θ and Γ_S are the non-dimensional transfer coefficients of heat and salt, respectively, determined from the near-wall physics. Based on a high-resolution LES study on the heat and salt transfer coefficients, which are described as functions of the friction velocity and thermal driving, these coefficients, Γ_θ and Γ_S, at the ice-shelf base were found to be 8 × 10⁻³ and 2.6 × 10⁻⁴ for basal melting, respectively (Vreugdenhil and Taylor, 2019). For the sea-ice formation at the sea surface, the same coefficients (5.8 × 10⁻³ and 2 × 10⁻⁴) used in a previous study on sea-ice formation in polynyas were used, because the thermal driving in our study was comparable to that in the previous study (Heorton et al., 2017). These melting and freezing effects were applied at the first grid from the ice-shelf base or the sea surface. In this study, the melting or freezing effects at the vertical side of the ice front were not included.

2.3 Simulation description

In this study, we conducted four LES simulations of the turbulence intensity within the IOBL to reveal the influence of turbulence on IOBL characteristics. Using the masking method, an idealized ice shelf with a depth of 280 m was described at the upper-left part of the simulation domain to indicate the ice-shelf depth near the front of the NIS (Briscolini and Santangelo, 1989; Stevens et al., 2017). Left (right) boundary was set to inflow (outflow) boundary condition. Moreover, initial conditions of velocity, temperature, and salinity were the same as inflow boundary conditions. For zonal velocity, potential temperature, and salinity, the inflow boundary condition was set to Dirichlet boundary conditions with temporally constant vertical profiles. The outlet boundary condition was determined to match the radiation boundary condition (similar to extrapolation), which prevented numerical errors without rapid changes in the velocity and scalar properties. This radiation boundary condition at the outlet boundary allowed wave-like motions within the domain to pass through the boundary with only a small reflection. The cyclic boundary condition was applied to lateral boundaries, whereas the Neumann boundary condition for momentum was imposed on the top boundary condition. In other words, the wind affects the scalar (temperature and salinity) fluxes but not the momentum fluxes at the sea surface. Details on the simulation domain and the boundary conditions are presented in Fig. 1, which shows the target study region, observation points, and simulation domain with a schematic diagram for the oceanic flow alongside a sub-ice-shelf plume. In this study, we used the theoretical profiles of velocity for describing the outflowing of the sub-ice-shelf plume. Different turbulence intensities were described in four different vertical profiles of the zonal velocity and imposed as inflow conditions at the boundary of the domain. Based on the power-law assumption of turbulent boundary layer flow ($U=U_{\mathrm{t}}\times (z/z_{\mathrm{0}}{)}^{(\mathrm{1}/n)}$), different velocity profiles (inset profiles in Fig. 2) were composed using different power-law indices, n=3 (weak turbulence) and 4, 5, and 7 (strong turbulence), to simulate the turbulence intensity within the IOBL (Irwin, 1979; Kikumoto et al., 2017) (height z represents the vertical distance from the ice-shelf base). A surface roughness (z₀) was 0.005 m (Gwyther et al., 2016). The freestream (geostrophic) velocity U_t was set as 0.06 m s⁻¹ at 572 m depth, based on the in situ observations near the ice front. To explore the effects of the sub-ice-shelf plume on ocean dynamics without wind stress, inflow boundary condition for velocity was set to zero above 280 m depth, and zero surface stress at the sea surface was used as the top boundary condition, whereas that below the ice-shelf base (280 m) was used for the sub-ice-shelf plume. At the start of the simulation, the initial profiles were imposed on the whole domain. The vertical profiles of the inflow boundary conditions for potential temperature and salinity were determined from the 24 CTD observations (Fig. 5). The potential temperature at inflow boundary condition above 280 m depth and at a depth from 280 to 570 m was set equal to the temperature at the sea surface (−1.9 ^∘C) and to the average temperature (−2.06 ^∘C), respectively. Below 570 m, the inflow profile of potential temperature was set to averaged potential temperature in CTD observation. The averaged salinity values (obtained from CTD observation) were adopted as the inflow boundary condition of salinity. The simulation dimensions were 3456 m × 3456 m × 864 m in the x, y, and z directions, respectively. For the simulations, a grid of 288 × 288 × 144 cells was used with a 12 m horizontal grid and a 6 m vertical grid. The grid size and surface roughness were adopted based on a numerical study on high shear plume flow (Gao et al., 2019). We performed a grid sensitivity study of the optimal grid size for higher accuracy and less computational costs (Fig. S1 in the Supplement) and found that a moderate grid resolution of 288 × 288 × 144 was suitable for resolving the turbulence with the parameterization of melting and freezing effects. The total simulation time required to reach the quasi-steady state was 96 h. A random generator for small velocity perturbation was applied at depths from 30 to 800 m to quickly spin up the turbulence. In the model simulation, we configured two physical domains, i.e., the sub-ice-shelf and open-ocean regions. In the open ocean, the simulated ocean velocity, potential temperature, and salinity results were validated using the CTD and LADCP observational data. Subsequently, the flow characteristics of the IOBL flow beneath the ice shelf were investigated in detail.

Figure 2 Time series of the friction velocities in all the four cases. The total time (96 h) was normalized by each large-eddy turnover time, t^∗, which was calculated from the overturning eddy scale (ice-shelf–ocean boundary layer (IOBL) scale) divided by the friction velocity. The inset figure shows four different theoretical profiles of the velocity.

3 Results

3.1 Quasi-steady ocean environment near the ice front

To confirm that the IOBL and oceanic flows approach a quasi-steady state, we plotted the time series of the friction velocities at the first grid below ice-shelf base (Fig. 2). The total simulation time (96 h) was normalized by the large-eddy turnover time (t^∗), which was calculated as the scale of overturning large eddy within the IOBL (IOBL depth) divided by the friction velocity. The friction velocities in the four LES cases showed significant fluctuations during the whole period, indicating that large-scale eddies exist beneath the ice shelf. The friction velocity fluctuations due to large-scale eddies show a repetitive pattern. The averaged friction velocities after 14 t^∗ for n=3, 4, 5, and 7 were 0.00169, 0.00208, 0.00255, and 0.00283 m s⁻¹, respectively (Ramudu et al., 2018). This difference is due to the different momentum entrainment by the turbulence intensity within the IOBL. We concluded that these IOBL flows approached a quasi-steady state after 14 t^∗ because averaged friction velocities for the four cases changed slowly enough. Table 2 presents these friction velocities and the large-eddy turnover times for the four cases. We calculated the time-averaged results within the last 3 t^∗ period for the later analysis to capture the averaged features of the flow without temporal variance.

Table 2 Main values of four different cases ($R{i}_{\mathrm{f}}=\frac{\frac{g}{{\mathit{\rho }}_{\mathit{\theta },\mathrm{0}}}\left(\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathit{\theta }}^{\prime }{w}^{\prime }}\right)}{-\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\left(\frac{\partial U}{\partial z}\right)}$).

Figure 3 illustrates the vertical sections (y=1728 m, domain center) of the zonal velocity, potential temperature, and salinity, which are time-averaged in the final 3 t^∗ period to examine the spatial distributions of the flow structures and variables in two end-member cases (n=3, n=7) for turbulence intensity. In the open-ocean region, the velocities for weak- and strong-turbulence cases (upper row of Fig. 3) exhibited similar patterns for the two overturning cells in the upper ocean region (0–280 m depth). Since we did not impose the wind effect at the top boundary, we can conclude that the development of the outer overturning cell is mainly induced by the downwelling (negative buoyancy flux) of the locally concentrated salinity as well as the shear stress induced by the momentum difference between the upper region and sub-ice-shelf plume. Moreover, the development of the inner overturning cell is mainly due to the upwelling of the buoyant water and the downwelling of the salt flux. Thus, the inner overturning cell is stretched in the vertical direction, whereas the outer overturning cell is stretched in the horizontal direction. This is because the driving forces of the inner cell and outer cell were the buoyancy force and shear force, respectively. The downwelling convection which the two overturning cells share drives the sub-ice-shelf plume down, moving the isopycnal line (280 m depth at initial state) near the ice shelf to a depth of approximately 350 m. Positive buoyancy at the left side (near the ice front) of the inner circulating cell originated from the meltwater created at the ice-shelf base near the ice front. In this study, we refer to this water as the positively buoyant ice-shelf water (PISW). In contrast to the sub-ice-shelf plume, which has a neutral buoyancy near the ice front, PISW has a strong buoyancy. Consequently, PISW is a major contributor to the formation of the inner overturning cell. Because the temperature of the sub-ice-shelf plume is higher than the local freezing temperature at a depth of 280 m, basal melting occurs (Fig. 7) and creates the PISW. At zonal distances from 1280 to 1600 m, this PISW mixes with the outer ocean and exhibits a temperature that is approximately 0.1 ^∘C lower than the surface freezing temperature, affecting the sea-ice formation near the ice front (middle row of Fig. 3). As shown in the lower row of Fig. 3, the upwelling of PISW causes an upward movement of the isopycnal line (identified by potential density). Except for the upper open ocean where overturning cells are dominant, the water column is stratified well below a depth of 350 m.

Figure 3 XZ cross-sectional contours (y=1728 m, domain center) of the zonal velocity, potential temperature, and salinity in the n=3 (left) and n=7 (right) cases. In these contours, the zonal direction is perpendicular to the ice-shelf front. These results are time-averaged in the last 3 t^∗ period. Upper row: zonal velocity; middle row: potential temperature; lower row: salinity. In the lower row, the isopycnal lines are identified by the potential density (kg m⁻³) with 0.01 intervals.

A noticeable difference between the two cases is observed near the ice front and beneath the ice shelf. At depths from 280 to 320 m (IOBL region), relatively high zonal velocity beneath the ice shelf is observed in the strong-turbulence case (n=7). After passing the ice shelf, this relatively high speed current flows in a direction perpendicular to the ice front. In the weak-turbulence case (n=3), a relatively low velocity is observed beneath the ice shelf, and the current rises after passing the ice-shelf front. These momentum differences in the two cases mainly affect the magnitude and scale of the circulating cells near the ice front. Another noticeable difference between the two cases is the different temperature and salinity in the upper ocean region. In the strong-turbulence case, the increased friction velocity affects the basal melting rate and a large amount of PISW, leading to a low temperature and salinity in the upper ocean region.

To examine the sea-ice distribution and PISW effect on sea-ice formation, we illustrate the horizontal distributions of the 3 t^∗ time-averaged freezing rate (sea-ice formation) at the sea surface, as shown in Fig. 4. Because the atmospheric forcing, friction velocity, interfacial temperature, and interfacial salinity at the sea surface are the same in all the four cases, the difference in the freezing rate originates only from the ocean conditions. Although the amount of PISW is small in the weak-turbulence case, most of its PISW is rising along the ice front owing to a low flow momentum within the IOBL. Otherwise, some of the PISW in the strong-turbulence case is also rising along the ice front, but most of it is advected to the upper mixed layer in the open-ocean region. As a result, different patterns of the freezing rate are observed in the two cases. In the weak-turbulence case, most of the freezing (sea-ice formation) is concentrated at the frontal region of the ice shelf, indicating a maximum freezing rate of 8.9 m yr⁻¹. Although a maximum freezing rate of 7.16 m yr⁻¹ is observed right in front of the ice shelf in the strong-turbulence case, the spatial-averaged freezing rate is higher than that in the weak-turbulence case. An interesting point in the strong-turbulence case is the heterogeneous pattern of the freezing rate in the meridional direction. This feature is strongly related to the PISW layer formed by the PISW upwelling (Fig. S2). In the weak-turbulence case, the PISW upwelling occurs along the ice-front edge, forming a strong and narrow PISW layer near the ice front with a strengthened inner overturning cell. However, heterogeneous patterns of the freezing rate are observed in the strong-turbulence case, because the PISW layer near the ice front is wide zonally with a weakened inner overturning cell, permitting the larger baroclinic disturbance caused by sloped isopycnals. This heterogeneous pattern of the freezing rate is comparable to the disturbance scale (2066 m), as identified from the Rossby radius of deformation, which indicates that the length scale the rotation effect is dominant. This scale is obtained from the depth-averaged buoyancy frequency and depth between the sea surface and IOBL bottom.

Figure 4 XY horizontal distribution of the 3 t^∗ time-averaged freezing rate (m yr⁻¹) at the sea surface. (a) Freezing rate in the case of n=3, weak turbulence. (b) Freezing rate in the case of n=7, strong turbulence. (c) Zonal–spatial distribution of the freezing rate in the four different turbulence cases. These values are averaged along the meridional direction.

3.2 Validation of simulation results

The CTD and LADCP observation data (grey lines in Fig. 5) shows the existence of sub-ice-shelf plumes, which exhibit lower temperatures than the surface freezing temperature. Above a depth of approximately 350 m, a well-mixed feature of potential temperature and salinity is observed, suggesting a strong mixing in this region. Two interesting observations that are difficult to explain were obtained. One of them is the existence of relatively low temperature freshwater (−1.96 ^∘C, 34.65 psu, practical salinity unit) at a depth of 100 m (Fig. 5b and c). Because of the sea-ice formation (latent heat and salt flux) at the polynya of the NIS region in the late summer season, we suggest that the low-temperature freshwater is produced from the ice-shelf base and not from the sea surface. This feature can be explained by the PISW upwelling process which can be observed in LES results. It is observed that the PISW is advected to the open-ocean region at 100 m depth in the potential temperature distribution at an instantaneous time of 21.8 t^∗ (not shown). Another is the negative velocity (toward the ice front) above a depth of 300 m. It is a coherent feature because all the 23 LADCP observations indicated a similar feature. Because the direction of the katabatic winds is positive (away from ice front), it is not the cause of the negative velocity. We hypothesize that there is a large-scale overturning cell (counterclockwise) caused by the shear force of the sub-ice-shelf plume and downwelling by the salinity flux at the sea surface.

Figure 5 Vertical profiles of the velocity, potential temperature, and salinity obtained from the CTD and LADCP observational data (solid grey lines), initial profiles (black dashed line), and the results of the large-eddy simulation (LES). (a) Zonal velocity; (b) Potential temperature; (c) Salinity.

To validate our LES results and examine our hypothesis for the large-scale overturning cell, we plotted the vertical profiles of the stream-wise zonal velocity, potential temperature, and salinity at a distance of 1 km away from the ice front and CTD and LADCP observations to compare the vertical distribution of the momentum and variables related to potential density (Fig. 5) quantitatively. In terms of the velocity, the LES results simulate the development of the outer overturning cell well, having vertical profiles similar to those observed in the LADCP data, although the velocity gradient of LES results is less sharp compared to that observed in the LADCP data. This mismatch between the LES results and the observations arises from the underestimated convection force (underestimated strength of overturning cells). The potential temperature and salinity profiles of these four LES results agree with the 24 CTD profiles, in terms of magnitude and depth of peak values. The difference in the potential temperature above a depth of 350 m also arises from the underestimated strength of the overturning cells. These LES experiments were initiated from an initial state with observation-constrained boundary conditions (grey lines in Fig. 5) and the effects of basal melting, sea-ice formation, and outflowing momentum of the sub-ice-shelf plume. Therefore, the underestimation of convection force comes from the setup of boundary conditions and from resolving these effects. However, we conclude that the LES results are consistent with the in situ observations for the oceanic environments, such as the development of the overturning cells and similar vertical structure of temperature and salinity.

Because the velocity gradient between the freestream velocity at 572 m and the velocity at the sea surface is similar in all the cases, the velocity profiles of the LES results at 1 km away from the ice front are similar. However, the different turbulent intensities affect the different momentum transfer into the IOBL, resulting in different melt rate. Due to a difference in the melt rate, the magnitude of the potential temperature in the upper mixed layer in all the cases is significantly different: in the upper mixed layer, the potential temperatures in the four cases are approximately −1.925, −1.93, −1.932, and −1.943 ^∘C. This difference is due to the amount of PISW created from the ice-shelf base and not due to the differences in the advection of PISW, because the total potential temperature in the strong-turbulence case is lower than that in the weak-turbulence case (Fig. 4). Similar features with the effect of PISW are also observed in the salinity profiles.

In the LES model, the filtered Navier–Stokes equation is solved with the modeled effect of small-scale eddies to reduce the computational costs. Therefore, the criteria for “small scale” are important; these criteria are determined by the grid size in the LES model. To evaluate and validate the grid size and the parameterization of small-scale eddies, it is necessary to confirm that the turbulence characteristics of the LES result are similar to the turbulence characteristics of the inertial subrange in which energy cascading occurs. We obtain the one-dimensional turbulence energy by integrating the inner product of the wavenumber and the two-point correlation calculated along the single line in the y direction. The one-dimensional turbulence energy spectra at a depth of 291 m (within the IOBL) in weak- and strong-turbulence cases are plotted in Fig. 6. Moreover, we examined different zonal locations (x=400, 800, 1200, 1600, 2000, and 2400 m) to observe a spatial transition of the IOBL flow. For the high wavenumbers which are comparable to the grid scale, the energy spectra slope of the LES results is slightly smaller than $-\mathrm{5}/\mathrm{3}$ slope of the Kolmogorov scaling in the inertial subrange. These are similar to the −1 slope of the Batchelor regime, indicating that the resolved turbulence has anisotropic characteristics in the high-Schmidt-number flows. For the weak-turbulence case, spatial transition of turbulence energy at the IOBL region is observed. As the direction of the flow is away from the inflow boundary, the turbulence energy is gradually increased beneath the ice shelf. In the open ocean, turbulence energy is similar. For the strong-turbulence case, the trend of energy spectra is similar to that for the weak-turbulence case, except at x=1200 m. The highest turbulence energy is observed near the ice front (x=1200 m), showing high momentum transfer near the ice-front edge.

Figure 6 One-dimensional turbulence energy spectra at a depth of 291 m at the PISW within the IOBL. (a) n=3 and (b) n=7. Different shapes and colors represent the values at different zonal distances: 400, 800, 1200, 1600, 2000, and 2400 m. These power spectra are obtained by y-direction (meridional direction) spatial assessment of time-averaged velocity at each x location. The $-\mathrm{5}/\mathrm{3}$ slope (Kolmogorov scaling) represents the regime of inertial subrange, whereas the −1 slope (Batchelor) represents viscous-convective range in the high-Schmidt-number flows.

3.3 IOBL characteristics in the sub-ice environment

The aforementioned analysis shows that the LES model adequately resolves the oceanic flow beneath the ice shelf with the thermohaline dynamics, such as IOBL dynamics, PISW upwelling, and convective downwelling by the salt flux at the sea surface. Next, we explore the flow characteristics of the IOBL beneath the ice shelf using the LES results. Since we assumed a flat base of the ice shelf in this study, the buoyant force of PISW does not accelerate the PISW at the ice-shelf base. The driving forces within the IOBL flow are the shear forces caused by the momentum of the inflow through the boundary and the stratification force (stabilizing force) caused by melting. In this section, we comprehensively analyze the oceanic flow characteristics to reveal the relation between the flow physics and melting patterns within the IOBL beneath the ice shelf.

Figure 7 illustrates the horizontal distributions of the 3 t^∗ time-averaged meting rate at the ice-shelf base. Evidently, different magnitudes of the melt rates are observed in the weak- and strong-turbulence cases. Except the magnitude of the melt rate, the overall trend in the two different turbulence cases is similar. In the region close to the inflow boundary and ice front, the melt rate is highest. Because the inflow boundary conditions are theoretical profiles, the transition for a fully developed (spatially quasi-homogeneous) IOBL flow is needed. As shown in Fig. S3 in the Supplement, turbulence intensities in the weak- and strong-turbulence cases are fully developed after 312 and 336 m distances, respectively. Therefore, we exclude the region with non-developed turbulence in the analysis. The spatially averaged values of melt rate in the two different turbulence cases are 0.092 and 0.153 m yr⁻¹, respectively, showing a 66 % difference in the melt rate near the NIS front. The melt rates obtained in this study are significantly lower compared to those reported by Wray (2019) (0.45–0.95 m yr⁻¹) and estimated via the CryoSat-2 satellite observation during 2010–2018 (1 ± 0.6 m yr⁻¹) at the NIS ice-front region (Adusumilli et al., 2020). The difference in melt rates can be explained by a discrepancy between our thermal driving, 0.056 ^∘C (2.116–2.06 ^∘C), and the thermal driving, 0.14 ^∘C (2.0–1.86 ^∘C), considered by the study of Wray (2019). In this study, only the effect of the sub-ice-shelf plume was considered, but the observations in summer season are likely affected by both intrusion of relatively warm Antarctic surface water and the effect of the sub-ice-shelf plume, resulting in a difference in the thermal driving and melt rates. If the melt rate in this study is assumed to be 0.12 m yr⁻¹ (averaged value of 0.092 and 0.153), we can estimate that 12 %–25 % of the total basal melting near the NIS front is due to sub-ice-shelf plumes.

Figure 7 XY horizontal distribution of the 3 t^∗ time-averaged melt rate (m yr⁻¹) at the ice-shelf base (280 m). (a) Melting rate in the case of n=3, weak turbulence. (b) Melting rate in the case of n=7, strong turbulence. (c) Zonal–spatial distribution of the melt rate in four different turbulence cases. These values are averaged along the meridional direction.

To observe and quantify the velocity structure within the IOBL in the strong-turbulence case and the Ekman layer development, we plot depth profiles (meridional-averaged) of the velocities, turbulence intensity, and horizontal momentum flux beneath the ice shelf (Fig. 8). At 280–400 m depths, the zonal velocities exhibits different structures for all the four cases. The strong-turbulence case displays the smallest mean shear gradient but the largest turbulence intensity, whereas the weak-turbulence case has opposing features. Since the turbulence kinetic energy production is proportional to the mean shear gradient and turbulent shear stress (Pope, 2000), the turbulent shear stress is the highest in the strong-turbulence case, indicating that the turbulent shear stress produces a large portion of turbulent kinetic energy production. Due to high turbulent shear, a strong frictional Ekman layer with negative meridional velocity (flows to the right of geostrophic flow) develops within the IOBL. Moreover, the frictional Ekman layer depths for the weak- and strong-turbulence cases are 11 and 17 m, respectively. These depths are comparable to the depths (11.9 and 19.6 m, respectively) estimated based on the friction velocity and Coriolis parameter (Coleman et al., 1990).

Figure 8 Vertical profiles of (a) u and v mean velocities and (b) turbulence intensity. The inset in panel (a) shows the vertical profile of the horizontal momentum flux for the Ekman layer formation beneath the ice shelf.

Figure 9 Vertical profiles of the momentum fluxes and heat fluxes. The fluxes were characterized using resolved flux and subgrid-scale (SGS) flux. (a) Momentum fluxes (resolved: ${\mathit{\rho }}_{\mathrm{0}}\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$; SGS: ${\mathit{\rho }}_{\mathrm{0}}{K}_{\mathrm{m}}\stackrel{\mathrm{‾}}{\partial u}/\partial z$). (b) Heat fluxes (resolved: ${\mathit{\rho }}_{\mathrm{0}}{c}_{\mathrm{s}}\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}^{\prime }}$; SGS: ${\mathit{\rho }}_{\mathrm{0}}{c}_{\mathrm{s}}{K}_{\mathrm{h}}\stackrel{\mathrm{‾}}{\partial \mathit{\theta }}/\partial z$, where ρ₀ (1024 kg m⁻³) is the reference density of seawater and c_s (4.02 × 10³ J kg⁻¹ K⁻¹) is the specific heat capacity of seawater (Sharqawy et al., 2010)).

Figure 9 shows the vertical profiles of the vertical fluxes of momentum and heat beneath the ice shelf. As shown in Fig. S4, which depicts the vertical buoyancy flux profiles, the PISWs in both the weak- and strong-turbulence cases have a stabilizing effect because of their positive buoyant characteristics. Moreover, the PISW in the weak-turbulence case exhibits a larger buoyant force than that in the strong-turbulence case. Combining the regions, where the meltwater and its stabilizing effect dominate, with that where the turbulence-induced heat entrainment is vigorous (5 % of the maximum heat flux), the IOBL bottom is determined, since the IOBL is analogous to the atmospheric boundary layer (Derbyshire, 1990). That is, the IOBL region was defined as the region where the thermodynamic changes induced by the basal melting at the ice-shelf base are dominant. In the strong-turbulence case, the vertical momentum flux is negative, and its maximum is located within the IOBL (IOBL bottom: 319 m depth). This implies that the momentum entrainment from the sub-ice-shelf plume to IOBL is effective, resulting in a large heat entrainment. However, the depth of the maximum negative flux in the weak-turbulence case is located at 347 m, i.e., slightly away from the IOBL. This difference causes the difference in the heat flux magnitude at the PISW bottom. Because the IOBL flow is quasi-steady, a similar temperature within the IOBL is maintained under heat entrainment from the sub-ice-shelf plume and cooling effect of the PISW advection. The positive heat flux at 280–319 m depths and negative heat flux at 320–400 m depths are caused by the large-scale turbulence convection, showing the occurrence of heat entrainment from the sub-ice-shelf plume and cooling effect of the PISW advection. The integrated area of heat flux represents the heat entrainment for the basal melting and PISW formation. The maximum positive heat fluxes for the weak- and strong-turbulence cases are 138 and 213 W m⁻², respectively, with a 54 % difference. This difference is comparable with that (66 %) in the melt rate near the ice front, confirming that the basal melt rate is proportional to the amount of heat flux and entrainment, because of the flow advection penetrating the stratified IOBL.

4 Discussions

4.1 Ocean environment near the NIS front

In this study, we simulated the neutrally buoyant, sub-ice-shelf plume and the oceanic environment affected by this plume. The main findings and conclusions are summarized in Fig. 10. First, we observe the development of two overturning cells near the ice front. The inner overturning cell is caused by the upwelling of the PISW along the ice-front slope and downwelling at the concentrated salinity flux (local salt maximum) induced by the sea-ice formation. The outer overturning cell is caused by outflowing momentum of the sub-ice-shelf plume. Second, we identify the role of turbulence within the IOBL. Notably, a higher turbulence intensity results in a large amount of ice melting near the ice front. A high turbulence intensity causes a large momentum transfer, resulting in an increased melting and relatively high speed currents within the IOBL. The large momentum causes the IOBL current (positive zonal velocity) to flow perpendicular to the ice front after it passes the ice shelf.

The horizontal scales of the overturning cells differ slightly for the weak- and strong-turbulence cases. The horizontal scale (492 and 408 m for the cases with weak and strong turbulence, respectively) of the inner overturning cell is calculated based on the locations of the local salt maximum (zonal distance = 1776 and 1692 m for the case with weak and strong turbulence, respectively) at the sea surface. Because we observed the negative velocity after local salt maximum at the sea surface, positive velocity of the sub-ice-shelf plume, and downwelling of the salinity flux, we can conclude that an outer overturning cell exists. However, we cannot observe the exact horizontal scale of the outer overturning cell owing to the limitation of the domain scale in this study. Because the flow with the sub-ice-shelf plume beneath the ice shelf is a backward-facing step flow (reattachment flow with geometry), we can estimate the horizontal scale of the outer overturning cell based on the reattachment length (Rygg et al., 2011). For the oceanic flow at a high Reynolds number (2 × 10⁶), we observe the development of the outer cell (5.3–8.4 h (h is geometry height)) and inner cell (1.3–1.5 h). As the ice-shelf thickness is 280 m, we suggest a 364–420 m inner overturning cell and 1484–2352 m outer overturning cell.

Figure 10 Final conclusions and schematic diagram representing the oceanographic picture near the NIS front and the various phenomena that occur within the IOBL, as resolved via LES.

PISW upwelling is important for the formation of modified Antarctic surface water and ocean circulation near the ice front. In a previous study on the freshwater wedge in the western Ross Sea (Malyarenko et al., 2018), the freshwater layer at sea surface was observed. They hypothesized that this freshwater mainly originates from the sea ice melting, ice-shelf water (ISW) outflow, and frontal ablation. The PISW upwelling and its accumulation observed in this study could evidentially support the “basal melting input” as a feasible wedge formation mechanism proposed by Malyarenko et al. (2019). If we consider the relatively warm Antarctic surface water with the inclusion of frontal ablation in a future study, we can investigate and quantify the contributions of the basal melting and frontal ablation in the freshwater input in the NIS region. In the study of meltwater outflow and its vertical structure by Garabato et al. (2017), they observed that the injection of the high-buoyancy, meltwater-rich glacially modified water triggers overturning via centrifugal instability near the ice front. For the vertical velocity of the PISW (Fig. S5), we observed a high positive velocity (0.04 m s⁻¹) near the ice front and a negative velocity (0.02 m s⁻¹) at the local salinity maximum. Given that approximately 0.02 m s⁻¹ of the vertical velocity is observed at 1 km away from the ice front by Garabato et al. (2017), our findings can support the ascent of meltwater outflow observed in this previously reported study. In contrast to centrifugal instability in the Pine Island Ice Shelf, we can demonstrate the symmetric instability in this study (Figs. S6 and S7). This difference is caused by the different directions of the current near the sea surface and the exclusion of the katabatic wind effect.

4.2 Limitations of the LES experiments

As shown in previous LES simulation of the ice–ocean boundary (Ramudu et al., 2018; Vreugdenhil and Taylor, 2019), the LES model is a powerful tool for resolving the flow and its three-dimensional structures with the change of mass and energy, with a high accuracy and low computational cost. Therefore, employing the numerical approach such as a regional ocean model and the LES model is one of the best solutions for resolving the three-dimensional structures in the oceanic flow. Under these circumstances, we demonstrated one method for investigating the IOBL physics and ocean dynamics by combining the numerical approach with observational data and theoretical profiles (power-law profiles with different turbulent intensity). In this study, we used the LES model with in situ boundary conditions to expand the one-dimensional observation profile in the open-ocean region to the three-dimensional flow-field in the open-ocean and sub-ice-shelf regions. Additionally, we set the interfacial values (θ_b, S_b) in the liquidus condition using ambient values (θ_f, S_f) obtained from the CTD observational results. In this study, to evaluate the applicability of our proposed method, we validated our results with observations, in terms of overturning cell characteristics, vertical structures of temperature and salinity. In a future work, this proposed methodology will be validated using in situ sub-ice-shelf observations made using advanced instruments such as hot water drilling, autonomous underwater vehicles, and glider, among others. Moreover, the boundary conditions for the sub-ice-shelf region, the effect of wind stress, and salt flux changes due to precipitation and evaporation will be investigated for simulating the IOBL and overturning cells more realistically.

In the Terra Nova Bay, the northeastward katabatic wind is dominant, and it drives the along-front current (Guest, 2021; Malyarenko et al., 2019). If we included the wind effect at the sea surface, the horizontal mixing may be enhanced, advecting the fresh meltwater of the PISW layer near the ice front to the open sea. In terms of overturning cells, the strength and horizontal scale of the outer overturning cell may decrease, as the wind stress reduces the shear stress between the sea surface and the sub-ice-shelf plume. The weakened outer cell, in turn, weakens the inner (secondary effect) cell. However, the strength and scale of the inner overturning cell may be similar to that obtained in the present study, because the wind stress at the sea surface is imposed at the upper part (positive velocity region) of the inner overturning cell. In summary, if the wind stress effect is included, similar scale of the inner circulating cell and a decreased scale of the outer circulating cell near the ice front can be expected.

In this study, we did not include the melting or freezing effects at the vertical side of the ice shelf. Because thermal driving at a depth of 146 m on the vertical side of the ice shelf was almost zero, the effects of melting and freezing occurred below and above 146 m depth, respectively. Because 146 m was almost near the center of the ice front, the total change in the temperature and salinity, because of the melting and freezing effects at the entire ice front, was not significant. However, this feature in the ice mass change at the ice front might be related to the shape of the ice-front edge.

In terms of the frazil ice processes, we only considered the latent heat release and salt flux produced by the sea-ice formation at the sea surface. However, suspended frazil ice, formed by adiabatic cooling, can occur in the region of PISW upwelling at the ice front and upwelling of the sub-ice-shelf plume (in this study, this upwelling occurred outside the study domain). The effects of the frazil ice dynamics (e.g., crystal growth rate, nucleation, and concentration) on the sea surface or ice-shelf water (PISW and sub-ice-shelf plume) should be investigated, because the change in plume characteristics as well as the temperature and salinity is strongly related to the frazil ice dynamics (Galton-Fenzi et al., 2012; Rees Jones and Wells, 2018). Because the inclusion of suspended frazil ice affects the increased plume velocity and decreased plume density, the upwelling of the PISW can be strengthened. This increases the strength of the inner overturning cell. Owing to a high plume velocity and decreased plume density, only some precipitation of frazil crystals within the PISW occurs at the vertical side of the ice front. However, if the ice-front geometry is modified through an increased frazil ice formation by the freshwater wedge formation, precipitation on the modified geometry might increase, exhibiting a nonlinear effect on the change in the ice-shelf geometry (Smedsrud and Jenkins, 2004).

In this study, we employed theoretical power-law profiles of the velocity, which exhibit different turbulence intensities, because observational data on the vertical structures beneath the ice shelf are rarely available. In the flux Richardson number in Table 2, the relationship between the stratification and turbulent mixing is not clear, even though we imposed different turbulence intensities. This implies that the negative feedback (enhanced buoyancy fluxes) of the PISW increases as the turbulence intensity and its entrainment increase. In this study, we performed a preliminary assessment of the oceanic structures in a sub-ice-shelf environment. The constant model coefficient (C_m) in the SGS model is applied in the IOBL flow near the wall and the flow within the sub-ice-shelf plume to reduce the computational costs. In future studies, corrections to the coefficients will be required, based on high-resolution studies of near-wall turbulence. The main findings of this study can be used for understanding the IOBL physics and meltwater behavior near the sub-ice-shelf cavity. If direct observational data for the IOBL flow structures and turbulence characteristics in the sub-ice-shelf environment can be obtained, then this study can be improved by comparing LES results to observations and by introducing corrections into the ambient values and transfer coefficients for potential temperature and salinity. Other important factors that were not included in this study should be addressed and investigated in future studies such as the effect of ice-shelf bathymetry (shape and slope) and surface roughness on the turbulence characteristics, temporal variability of the sub-ice water plume, wind stress effects, and consideration of the melting at the vertical side of the ice front. A better understanding of the relationship between aforementioned factors and the IOBL physics with basal melting will help in improving the parameterizations (e.g., vertical mixing within the IOBL and sea-ice formation and behavior) in the regional ocean model (e.g., ROMS and MITGCM).

5 Conclusions

We successfully simulated the IOBL flow with a sub-ice-shelf plume, the melting effect beneath the NIS, and the freezing effect at the sea surface using a LES model and boundary conditions based on in situ observations. Since the inflow profile beneath the ice shelf is difficult to observe in in situ observation, we assumed the theoretical power-law profile to describe the different turbulence intensities. The flow simulated for a period of 96 h reached a quasi-steady state after 14 t^∗ (69 and 55 h for the cases with weak and strong turbulence), exhibiting a large variance and repetitive pattern with temporal convergence of friction velocity. To validate the LES results for the four different turbulence intensities, the simulated zonal velocity, potential temperature, and salinity were compared with the 24 h CTD and LADCP observational data obtained near the NIS in the Terra Nova Bay of the western Ross Sea, during a 2016–2017 shipboard survey. The simulation results are consistent with the observations considering the development of two overturning cells and thermohaline properties.

The different turbulence intensities in the IOBL flow beneath the ice shelf resulted in significantly different melting features and flow dynamics. With increased friction velocity, the melt rate in the strong-turbulence case increased by 66 % compared to that in the weak-turbulence case, maintaining the stratification intensity (similar flux Richardson number). In the strong-turbulence case, distinct features such as higher basal melting (0.153 m yr⁻¹), weak upwelling of the PISW, and a larger sea-ice formation were observed, suggesting a relatively high speed current within the IOBL because of a highly turbulent mixing. Comparing with the observations, we estimated that 12 %–25 % of the total basal melting near NIS front is caused by the sub-ice-shelf plumes. We observed that the sub-ice-shelf plume, PISW upwelling, and downwelling of concentrated salinity flux compose two overturning cells near the ice front.