Flow aware parameterizations invigorate the simulated ocean circulation under the Pine Island ice shelf, West Antarctica

Warm, subsurface ocean waters that access ice shelves in the Amundsen Sea are likely to be a key driver of high meltrates and ice shelf thinning. Numerical models of the ocean circulation have been essential for gaining understanding of the mechanisms responsible for heat delivery and meltrate response, but a number of challenges remain for simulations that incorporate this region. Here, we develop a suite of numerical experiments to explore how sub ice shelf cavity circulation and meltrate patterns are impacted by parameterization schemes for (1) subgrid-scale ocean turbulence, and (2) ice-ocean interactions. To provide a realistic context, our experiments are developed to simulate the ocean circulation underneath the Pine Island ice shelf, and validated against mooring observations and satellite derived meltrate estimates. Each experiment is forced with data-informed open boundary conditions that bear the imprint of the gyre in Pine Island Bay. We find that even at a  ̃600 m grid resolution, flow aware ocean parameterizations for subgrid-scale momentum and tracer transfer are crucial for representing the circulation and meltrate pattern accurately. Our simulations show that enhanced meltwater diffusion near the ice-ocean interface intensifies near wall velocities via thermal wind, which subsequently increases meltrates near the grounding line. Incorporating a velocity dependent ice-ocean transfer coefficient together with a flow aware ocean turbulence parameterization therefore seems to be necessary for modelling the ocean circulation underneath ice shelves in the Amundsen Sea at this resolution.

to explore how sub ice shelf cavity circulation and meltrate patterns are impacted by pa- informed open boundary conditions that bear the imprint of the gyre in Pine Island Bay. 28 We find that even at a ∼600 m grid resolution, flow aware ocean parameterizations for lead to ice shelf melting and increased glacial mass loss. Simulating the ocean waters that 42 reach ice shelves remains challenging, however, because it is difficult to accurately rep-43 resent turbulence in the ocean and interactions at the ice-ocean interface. In this study, 44 we show results that can provide practical guidance for accurately capturing these pro-45 cesses in computer models. We focus on developing a simulation of the ocean circula-46 tion under the Pine Island ice shelf, which is fed by one of the fastest flowing glaciers 47 in Antarctica. We show that when the speed of ocean currents is used to determine the 48 rate of heat and salt exchanges due to turbulence, the resulting simulation resembles ob- Our goal in this study is to provide understanding for the various parameteriza-116 tion choices available to ocean-only models that simulate the circulation under ice shelves 117 in the Amundsen Sea. As such, we develop a number of numerical experiments to test 118 how these parameterizations impact the cavity circulation under the Pine Island ice shelf. 119 The unique configuration for each experiment is shown in Table 1. Here we outline the 120 study region and general model configuration that is applicable to all experiments. 121 The computational domain includes the cavity underneath the Pine Island ice shelf. 122 It extends westward to 102.75 • W, and northward to approximately 74.46 • S, see Figure   123 1. Temperature, salinity, and zonal velocity is specified at the western open boundary,  free of sea ice during the simulated time period (Scambos et al., 1996), see section 2.4. 135 We specify the bathymetry and ice topography by regridding output from BedMachine  , where the minimum cell size is 2 m. 144 We remove ice from grid cells where the regridded ice topography is only < 0.2 m, 145 such that these grid cells are ice-free. We remove ice from these areas because the com-146 puted heat fluxes in these areas is unreasonably high, due to a division by the ice thick-147 ness. We note that the cutoff chosen here (0.2 m) is arbitrary, and we found values less  We approximate an initial condition for the model spinup by "extruding" the tem-      190 where C d = 1.5×10 −3 is the drag coefficient at the ice-ocean interface, and U M is the 191 near wall velocity: 193 Here (·) BL denotes a vertical volumetric average over a boundary layer that is one grid   Table 1. In these ex-213 periments the Laplacian and biharmonic viscosities are chosen to be a fraction of: based on the CFL criterion for numerical stability (Griffies & Hallberg, 2000), where L 216 is the local grid scale: (1)

218
With a nominal grid spacing such that L 600 m across the domain, and ∆t = 150 s 219 for these four experiments, the Laplacian (biharmonic) viscosity is roughly 120 m 2 /s (540,000 m 4 /s) 220 for base and constIO, and 18 m 2 /s (81,000 m 4 /s) for smallVisc. We note that the vis-221 cosity values for the first two experiments appear to be high, but are necessary for nu-222 merical stability in constIO. We therefore use the same values in base for comparison. 223 We specify only a Laplacian diffusivity for horizontal tracer transport, which is taken to the first baroclinic Rossby radius of deformation given by Chelton et al. (1998): Near the grounding line this ratio is below 2, such that the effect of the largest eddies 236 and baroclinic instabilities are only partially resolved (Hallberg, 2013). On the other hand, 237 farther away from the grounding line the resolution is well above the deformation radius.

238
Therefore, even at this sub-kilometer resolution, the model is in a gray zone, motivat- Leith (Table 1) are chosen for numerical stability. In these simulations, it is unclear how 248 to specify the diffusivity field, and we therefore tested the effect of various formulations 249 for the diffusivity tensor and intensity κ h . With a diffusivity tensor acting aligned with 250 the grid, we tested κ h = 1.0 m 2 /s, which made no discernible difference to κ h = 0.1 m 2 /s 251 in Leith. We additionally tested the effect of rotating the diffusion tensor along isopy-252 cnals as in Redi (1982), and found that this had a negligible effect on the resulting sim-253 ulation as well. of the GM scheme (Griffies, 1998), such that κ ρ = κ GM = ν h . The resulting diffusion 262 tensor is: where S x = −∂ x σ/∂ z σ and S y = −∂ y σ/∂ z σ are the isoneutral slopes and σ is the lo-265 cally referenced potential density. While this formulation implies a small vertical diffu-266 sivity, we add an additional background value of κ r = 10 −4 m 2 /s for numerical stabil-   Finally, we note that we also tested the effect of using the parameterization pre-  Secondly, we use the mooring data within the computational domain to validate 303 (or invalidate) the equilibrium state of the numerical experiments described previously. 304 For this task, we must choose a subset of the mooring data that is consistent with the 305 data that is used to obtain the open boundary conditions. Therefore, we select data taken  , 1996). This is consistent with our modelling assumption that sea ice is excluded. 311 We compute the temporal mean and standard deviation of the mooring data at each 312 instrument location during the time periods outlined above to obtain a representative 313 state we can compare our models to. At most depth levels the temporal standard de-314 viation, σ, is small, so we prescribe the minimum values: which provides a means of representation error, i.e. error due to misrepresentation of point 317 data within the model grid cells. Using these minimum values also accounts for poten-318 tial conflicts between different observed values that correspond to the same grid cell, or 319 nearby neighbors. More details related to raw mooring data processing, including con-320 siderations involved with computing potential temperature from in situ temperature, are 321 provided in Appendix A.    and salinity (green). Lower numbers imply a better fit, and the base-10 logarithm is shown to emphasize that experiments with a value less than zero imply that the data misfit is lower than the standard deviation.  parameterization. We therefore present the barotropic streamfunction underneath the 388 ice shelf to give a summarized view of the circulation in each case, Figure 6. 389 We note at the outset of this discussion that no model represents the broad pat-   The left column (Figure 7(a,d,g)) shows the difference in the total horizontal diffusive The right column, Figure 7(c,f,i), shows the velocity difference between the two exper- Here v ⊥ is the velocity normal to the section indi-

586
We use these values to account for measurement error and, more importantly, represen-587 tation error, accounting for spatiotemporally localized features that we cannot or do not 588 want to infer during the optimal interpolation. The potential temperature data show spu-589 rious jumps, see for example in Figure 3(e) at ∼200 m depth. These temperature fluc-polation cannot succsefully capture, and we do not wish to represent in our equilibrium-592 state model. As such, we choose a fairly large uncertainty to cover these cases. The salin-593 ity data shows no such jumps, so it seems reasonable to provide a relatively small un-594 certainty. Finally, we note that we only use data with the highest quality control flag, 595 but that this did not remove any data that we considered using. 596 Next, we describe the steps we took to prepare the mooring data for our study. We 597 first bin average the temporal data to hourly time stamps. All data from a single instru-598 ment are assumed to be at a single depth level. This assumption ignores temporal depth 599 variability, which we find to be reasonable because the amplitude of variability is well 600 below the vertical resolution of our grid (20 m).

601
Some moorings do not have salinity data, so in these cases we represent salinity at 602 these locations with data from the nearest CTD, which is <1 km away. In such instances, 603 we double the observational uncertainty of the salinity estimate, noting that this makes 604 it consistent with the CTD data described above: σ CT D,S = 2σ M,S . With in situ tem-605 perature and salinity at each mooring depth, we convert in situ temperature to poten-606 tial temperature using PyGSW (Campbell, 2012).

607
In the case of the PIG S mooring during 2014, data at some depth levels are in-  Here f : R Nm m → d ∈ R N d is simply a linear interpolation operator, mapping the parameter fields to the location of available data.
As a matter of computational convenience we make the following assumptions. First, we assume that each parameter field is independent from one another, allowing us to solve three OI problems for temperature, salinity, and velocity separately. Second, we assume that the observational and prior uncertainties can be described by Gaussian statistics.
We further assume that the observations are independent, such that Γ Obs = diag{σ 2 Observational uncertainties (standard deviations) are assumed to be 0.5 • C for potential temperature and 0.05 g/kg for salinity as they are not provided, see Appendix A in the main text for details. The LADCP velocity data is provided with uncertainty estimates, which we use.
We specify the prior covariance as Matérn class due to the link between Matérn class Gaussian fields and the solution of the elliptic stochastic partial differential equation October 19, 2021, 9:15pm where W(x) is a standard white noise process. We employ the empirical relationship The last ingredient is the initial guess for the OI problem, m 0 . Simple inspection of the temperature and salinity data shows that these fields have mostly vertical structure, with slight variations in the depth of thermocline and halocline due to their horizontal location.
Therefore, we specify θ 0 and S 0 as vertical profiles based on polynomial regressions of the data. We note that using this has similar results to specifying θ 0 = 0 • C and S 0 = 34.36 g/kg, but the former provides a better fit to the observations. The spatial structure of the velocity data is less obvious a priori and we therefore specify u 0 = 0 m/s. Given these assumptions and specifications, the minimization problem in equation (1) is linear and we can write the solution to each independent OI problem as: Here, potential temperature is shown as an example, and a similar solution is obtained for salinity and velocity. Before these results can be used directly as forcing for the ocean model, the spatial integral is removed from the zonal velocity: Removing the spatial mean ensures that we do not add or remove mass from the domain, and there is no artificial sea level rise during the spinup to reach equilibrium. In practice, October 19, 2021, 9:15pm X -4 : this corresponds to removing a small average velocity: 0.00943 m/s. The resulting fields are shown in comparison to the observational data in Figure S1.