Representing grounding line migration in synchronous coupling between a marine ice sheet model and a z-coordinate ocean model

Abstract Synchronous coupling is developed between an ice sheet model and a z-coordinate ocean model (the MITgcm). A previously-developed scheme to allow continuous vertical movement of the ice-ocean interface of a floating ice shelf (“vertical coupling”) is built upon to allow continuous movement of the grounding line, or point of floatation of the ice sheet (“horizontal coupling”). Horizontal coupling is implemented through the maintenance of a thin layer of ocean ( ∼ 1 m) under grounded ice, which is inflated into the real ocean as the ice ungrounds. This is accomplished through a modification of the ocean model’s nonlinear free surface evolution in a manner akin to a hydrological model in the presence of steep bathymetry. The coupled model is applied to a number of idealized geometries and shown to successfully represent ocean-forced marine ice sheet retreat while maintaining a continuous ocean circulation.


Introduction
A number of important physical processes in coastal oceanography involve the horizontal influx of water into regions that were previously "dry", as well as the complete removal of water from other regions that were at some point "wet". In estuarine regions, submerged boundaries change with the tidal cycle, and numerical codes which attempt to model important biological and geomorphological processes must capture this "wetting and drying" accurately (Hardy et al., 2000;de Brye et al., 2010). Flood models of storm surge must properly capture the advance and retreat of the flooding front (van't Hof and Vollebregt, 2005). There is an extensive numerical literature which deals with the problem of encroaching and retreating coastal flows (Medeiros and Hagen, 2013).
There is, however, an important coastal-oceanographic process not often discussed in the computational literature surrounding wetting and drying problems. In Antarctica (and to a lesser extent Greenland, and likely elsewhere during past glaciations), the ice sheet is marine terminating: it extends into the ocean in the form of large floating ice shelves. Due to the relatively small ice/ocean density differential, this occurs at depths ∼ 500-1000 m below sea level. The location where the ice sheet goes afloat, called the grounding line, is a topic of much discussion in the literature surrounding ice sheet dynamics. This is due to it being a sharp transition between two very different regimes of ice flow (Vieli and Payne, 2003;Pattyn et al., 2006;Schoof and Hewitt, 2013). From an ice dynamics perspective, determining the floatation point is equivalent to a viscous contact problem, one which requires very sophisticated numerical schemes (Schoof, 2011), although certain approximations make the problem more tractable (Goldberg et al., 2009;Cornford et al., 2013).
The focus of this paper, however, is not on the glaciological dynamics of grounding line migration, but rather on the coupling between the ice and the ocean underneath. The ocean circulates in a cavity bounded above and below by the ice shelf and bedrock, respectively, and on the landward side by the grounding line, where the water column depth pinches off. The circulation within the cavity is influenced by density variation, and by the topography of the ice shelf and sea bed (MacAyeal, 1984). From the ocean's perspective, the ice sheet/ ice shelf presents itself as variable surface pressure, and when the surface pressure favors the flooding of previously "dry" domain, the ocean will do so. Aside from the spatially varying surface pressure, the problem is analogous to run-up on a sloped beach.
This wetting/drying problem is quite an important one in the context of Antarctic Ice Sheet contributions to sea level rise. Melting and thinning of floating ice does not contribute significantly to sea level change (Jenkins and Holland, 2007), but retreat of the grounding line towards the interior represents loss of grounded ice, which does. Moreover, if the depth of the bed deepens inland (as is the case around much of the Antarctic coastline, Fretwell et al., 2013), retreat of the grounding line can lead to an increase in ice sheet thinning rates upstream, potentially leading to a positive feedback effect (Weertman, 1974;Schoof, 2007b;Vaughan and Arthern, 2007).
Yet most ocean models, including those adapted to study ice shelfocean interactions (e.g., Holland and Jenkins, 2001;Little et al., 2008;Walker and Holland, 2007;Gwyther et al., 2014), have not implemented such wetting and drying. There are two important differences between the ice sheet wetting/drying problem and the "standard" coastal wetting/drying problem. Firstly, while the latter can be addressed through shallow-water equations, the former must be addressed by a three-dimensional model representing baroclinic motions in the ocean, and therefore must represent active tracer (heat and salt) evolution, as well as an evolving upper boundary (the base of the ice shelf). Hence the approaches used for wetting/drying problems in coastal oceanography cannot be straightforwardly transferred to ocean models.
Secondly, the wetting front advance and retreat associated with grounding line change is much slower than in flooding and storm surge problems, with observed retreat rates of up to ∼ 3 km/year (Rignot et al., 2014) but often much slower. The separation of time scales allows for quasi-static ("discontinuous") approaches, in which a sequence of ocean runs are carried out with fixed cavity geometries, with the geometry sequence arising from the evolution of the ice sheet model (Grosfeld and Sandhager, 2004;Goldberg et al., 2012;De Rydt and Gudmundsson, 2016). Such an approach makes the assumption that a given geometry and far-field oceanic conditions will lead to a unique circulation, which may not be the case (e.g. if the ocean cavity is still in a transient state when the ice adjusts). Additionally, the ocean model must be spun up with each coupled time step, making the approach unsuitable for regional or global ocean models. More recently efforts have been made to retain the ocean "state" during adjustments to the cavity geometry. Such "asynchronous" approaches are still not ideal: typical coupling frequencies of a month to a year may lead to significant depth changes in the ice-ocean interface upon geometry updates. The necessary infilling with predefined water properties can lead to violations of mass and tracer conservation, as well as significant nonphysical adjustments. In one instance, the latter has been addressed through imposing that barotropic velocities remain fixed between updates (Asay- Davis et al., 2016). Still, it seems clear that a synchronous approach, i.e. one in which the ocean geometry is adjusted on or close to the ocean time step, is a preferable approach to modelling ice sheetocean interactions on continental and global scales, particularly if the ocean is subject to forcing on fast time scales, such as changes in wind stress (Christianson et al., 2016) and episodic additions of fresh water (Smith et al., 2017).
In this paper, we modify the Massachusetts Institute of Technology general circulation model (MITgcm, Marshall et al., 1997) to allow for synchronous coupling of an ice sheet and ocean model. MITgcm is a general-purpose fluid solver for simulating process-level to global-scale ocean circulation that is usually configured in hydrostatic and Boussinesq approximations in vertical z-coordinates (as has been done in the present study). Components of the development have been completed previous to this studythe most important of which allows continuous thinning and thickening of a floating ice shelf, which we term "vertical coupling" . Here we focus primarily on a scheme to allow for both grounded and floating ice and a dynamic grounding line ("horizontal coupling"). In the following, we briefly discuss the vertical coupling scheme and demonstrate how it can be used to allow for grounding line migration over a flat bed topography with minimal code changes. We then discuss the difficulties involved with variable bed topography (which are common to all z-coordinate models), and present a strategy to overcome them, which involves combining the ocean model algorithm with a scheme akin to flow through a porous medium. Finally, we present results from the first three-dimensional synchronously coupled ice-ocean model of marine ice sheet retreat.

Flat topography: methodology
With a flat bed topography, we are able to cleanly simulate synchronous coupled grounding line migration by making use of recent novel developments within MITgcm to allow for "vertical coupling", defined above. The work builds on previous developments to allow thermodynamic ice shelf-ocean interactions within MITgcm (Losch, 2008;Dansereau et al., 2014) as well as the development of an ice sheet component within the modelling framework (Goldberg and Heimbach, 2013). Vertical coupling is described in detail in Jordan et al. (2017), but we briefly describe those components which are relevant to our study in order to provide context for our results and for our further developments of the model. Note that we refer below to the − z coordinate implementation, not the z* implementation .

Vertical coupling
In this subsection we give details of vertical coupling and the MITgcm glacial flow model, which are also described in Jordan et al. (2017). Readers familiar with this paper might skip to Section 2.3.
Vertical coupling within MITgcm hinges on the nonlinear free surface capabilities of the model . The free surface elevation η defined relative to a reference surface elevation = z d, which for the ice-free ocean is = d 0, but for a cell occupied by the ice shelf is generalised to the height of the ocean-ice shelf, = − z d, in Losch (2008). η is updated in each time step in a fully mass-, heat-, and salt-conserving fashion, and responds both to barotropic pressure gradients and to gradients in surface load p surfwhich is imposed as the weight per unit area of the ice shelf. As flexural stresses within the ice shelf are not presently considered, under the ice shelf, where g is gravitational acceleration, ρ i is ice density and H is ice thickness, which is updated at each time step in response to ice dynamics and basal melting or freezing. The ice model, rather than updating its velocity and thickness at the same time (as is common practice in ice sheet modelling), updates its thickness on the ocean time step. (Velocity updates, which are more costly, take place every 12 hbut this is acceptable as velocity change induced by thickness changes on this time scale is very small.) In this manner, the surface load can be updated smoothly without exposing the ocean to sudden, large changes in surface pressure. As the ice and ocean codes are both components of MITgcm, there is no issue passing ice thickness to the ocean code and melt rate to the ice code.
In the z-coordinate free surface implementation of MITgcm, the height of the top-level cell grows with η (in contrast to the z*-coordinate implementation, in which cells thicken and thin uniformly in a column). Without intervention this can either lead to poor representation of the ice-ocean boundary layer (Jenkins, 2016) as the ice thins, or to negative cell height as the ice thickens. To this end a "remeshing" algorithm has been implemented. Upon initialisation of MITgcm, model cells are flagged as being either ice or ocean. The remeshing process essentially allows cells to switch from ice to ocean, and vice versa, within a model run and without the need to reinitialise ice and ocean masks. Whilst the topmost ocean cell thickness in a given column evolves every time step, at predetermined intervals we check to see if it has grown above a "splitting threshold" or below a "merging threshold". In either event, a remeshing process is triggered, whereby cells that are too thick are split into two separate cells, and cells that are too thin are merged with the cell below (see Fig. 1 for visualisation). In both cases, salt, temperature, mass and momentum are locally conserved. In the present study, the splitting and merging thresholds are 1.3 and 0.29 times the baseline cell thickness in all simulations.

Ice sheet/shelf model
The ocean model is coupled to the MITgcm ice sheet flow model, developed within the MITgcm code and described briefly in the following. Glacial ice is modelled as a slowly flowing viscous fluid with a non-Newtonian power-law rheology (Cuffey and Paterson, 2010). Due to its viscous nature, forces are at all times in balance and its velocity field is determined by its geometry , thus reducing the Navier-Stokes equations to a Stokes flow problem. The MITgcm ice model uses an approximate force balance, in which the forces within the ice sheet are considered to be hydrostatic, and "membrane stresses" (internal stresses which transmit force horizontally) are considered depth-uniform . It has been shown that such an approximation is appropriate where motion is dominated by either vertical shearing or along-flow stretching and cross-flow shearing (Schoof and Hindmarsh, 2010;Goldberg, 2011), which encompasses most flow regimes in Antarctica. The equations solved for velocity (see Goldberg, 2011 for further detail) are Here H is vertical ice thickness, S is surface elevation, = U U V ( , ) is horizontal ice velocity, and ρ i is ice density. Overbars indicate vertical averaging while the b subscript indicates basal velocities. A linear frictional sliding law is prescribed by (2) and (3) for which β 2 is the sliding coefficient (the exponent of 2 indicating it is strictly positive). The above equations differ for floating ice (ice shelves) only in that = β 0. The rheology is determined by Glen's Law (Cuffey and Paterson, 2010), wherein = n 3 and A is a parameter governed by the fabric and temperature of the ice. In this study, β is uniform (under grounded ice) and A is uniform as well. All relevant parameter values can be found in Table 1 where = n n n ( , ) x y is the seaward normal to the calving front, ρ ref is the ocean reference density and z low is ocean bathymetry. At all other lateral boundaries, we impose no-slip conditions on the ice model.
The ice thickness evolves as a result of horizontal ice mass flux divergence and surface and basal mass balance, the last of which takes place only where ice is floating. The evolution equation is where m is ice shelf mass balance (positive where ice is melting) and a surface mass balance (set to zero in this study). The grounding line is dynamic and is determined by a simple floatation condition, a consequence of the assumption of hydrostasy: and where ice is floating, the surface and basal elevation (B) depends only on thickness and reference densities: Note that the criterion (8) is not consistent with the ocean model criterion, i.e. that implied by M min in (15), and the basal elevation B is not necessarily equal to the upper surface of the ocean. In our simulations, though, the two floatation criteria are not seen to disagree by more than a grid cell. It would be straightforward to modify the ice flow model to enforce these consistencies, though for the present study we have not done so.
Ice sheet models generally evolve thickness and velocities with the same time step, but to evolve both on the ocean time step would be prohibitive. Rather a "split" approach is taken. It is seen from Eq. (4) that the equations for velocity are nonlinear and must be solved iteratively. In uncoupled mode, the MITgcm ice model solves a sequence of linear elliptic partial differential equations to a required tolerance; the number of iterations required depends on the change in thickness over a time step. In our coupled model we take a single such iteration once to twice daily. In Jordan et al. (2017), it is shown this approach maintains convergence of the velocity solution. Meanwhile ice thickness is evolved on the ocean time step; the most recent velocities are used in Eq. (7).
For brevity we omit details on the thermodynamic coupling between ice and ocean. The basic equations are given in Holland and Jenkins (1999) and have been implemented in MITgcm by Losch (2008) and Dansereau et al. (2014), and modified by Jordan et al. (2017) to allow for real, rather than virtual, fresh water fluxes. Aside from those listed in Table 1, all parameters relevant to the melt are as in Dansereau et al. (2014).

Parameterised buttressing
In some of the results in this paper, a two-dimensional ( − x z) flowline ice model is considered, coupled to a corresponding two-dimensional ocean model (i.e., Coriolis forces are omitted) to enable long model runs that are computationally cheaper. This removes the effect of buttressing, whereby confining stresses acting on the ice shelf are transmitted back to the grounding line, slowing grounded ice motion (Thomas, 1979). Ocean melting feeds back on ice sheet flow by thinning the ice shelf, which reduces buttressing. The lack of buttressing is dealt with through a parameterisation, which makes a further approximation that the along-flow normal stresses in the ice shelf do not vary across the shelf (Dupont and Alley, 2005). Assuming flow is in the y- Here W is the half-width of the ice sheet/ice shelf, treated as constant. In the three-dimensional simulation we present later, the buttressing parameterisation is not needed.

Subglacial layer
The method of vertical coupling described above has been thoroughly tested in the case of ice shelves with no grounded component . However, alone it cannot be used to simulate grounding line migration. The difficulties in implementing such a capability are linked to the way in which the ocean surface is evolved. Here, the present algorithm for free surface evolution in MITgcm is briefly discussed (for further detail see Marshall et al., 1997 andCampin et al., 2004). In one horizontal dimension, the free surface equation can be written Here the n and + n 1 superscripts represent the current and new time steps (not to be confused with n in Eq. 4). FW represents net fresh water flux (combined precipitation, evaporation, run-off, ice-shelf melting and freezing), as a mass flux in kg/s. h is ocean column thickness, u is velocity. The overbar represents a vertical average, while the superscript indicates both the update of velocity from time n due to inertia, baroclinic pressure, viscous effects, etc. (referred to below as + G u n 1/2 ), as well as a correction due to momentum imparted by the barotropic pressure gradient over a time step: Combining the two yields the following equation for + η n 1 : This elliptic partial differential equation (which in two horizontal dimensions requires a large linear system solver, such as Conjugate Gradients), arises from the use of + η n 1 in the right hand side of (13) rather than η n . This implicit approach damps high frequency barotropic waves and computational modes (Dukowicz and Smith, 1994). It also allows for tractable ocean simulation by greatly increasing the allowable time step. On the other hand, it causes difficulties transporting fluid into ocean columns which were previously dry (column thickness zero): if the equation were solved over portions of the domain with zero thickness, there is no guarantee that negative thicknesses would not arise.
A compromise to maintain the efficiency of an implicit free surface is to maintain a thin ocean layer under all regions of the ice sheet which could potentially ungrounda so-called "thin-film" approach to wetting and drying (Medeiros and Hagen, 2013). The only complication arising is the fact that the weight of a grounded ice column can be far greater than the weight of the ocean column it replaces, driving a strong pressure gradient which will quickly deflate such a layer.
To address this we modify the definition of p surf above. A "minimum" ice mass M min is defined as the minimum mass required for ice to be grounded: where z low is the same as that used for the ice flow model, ρ avg is the horizontally averaged density over all ocean-filled cells at a given grid level, and z sl is the ocean surface height averaged over the ice-free ocean. h mwct is the minimum water column thickness. The definition of p surf above is modified as follows: Accordingly, the load felt by the ocean is the actual weight of the ice column where ice is floating, i.e. where ρ i H < M eff , and is equal to M min D.N. Goldberg et al. Ocean Modelling 125 (2018) 45-60 under grounded ice, with a narrow transition between the two regimes around the grounding line, as shown for a simple example in Fig. 2. The result is that when ice is thin enough to float with an ocean column of at least h mwct , the full load of the ice is felt; above this, the load is reduced. Note that (17) is heuristically defined to achieve a certain effect, rather than arising from physical principles. We further impose that there is no melting where ρ i H ≥ M eff . The approach has already been used in another study . In all results presented, ε m1 is 0.01 and ε m2 is − 10 5 . The sensitivity of results to h mwct is tested later in the paper.

Flat topography: experiment
To illustrate the applicability of the subglacial water layer method described in Section 2.3, a simple two-dimensional ( − x z) coupled model with a flat topography is applied. The coupled ice-ocean model is initialised in a two-dimensional domain, implemented through free-slip ocean walls and setting the Coriolis parameter to zero. The domain is 320 km long with a horizontal resolution of 800 m. Bathymetry is −1000 m uniformly and vertical resolution is 10 m. The ice flows into the domain at = x 0 with a fixed volume flux of 4 × 10 5 m 2 /yr per unit width. An initial calving front is specified at x = 240 km; the ice front can retreat from this point by thinning to zero thickness, but cannot advance beyond it. The ocean model extends to x = 320 km, and is forced with a "sponge" layer extending over the last 20 grid cells of the domain. In the sponge layer, temperature and salinity are relaxed to the profile in Fig. 3(a), with a time constant varying from 2 days at the domain boundary to 10 days at the edge of the sponge layer. The target thickness of the subglacial water layer, h mwct , is 1 m. The ice sheet model is initialized without coupling to the ocean and run to a steady state; the coupled model is then initialised using this ice sheet geometry and with a time step of 200 s, and run until either a steady state, or complete ungrounding of the ice stream, is reached. As in Snow et al. (2017), velocities at the open-ocean boundary are adjusted in response to the spatially-averaged open-ocean surface elevation in order to maintain sea levels. Without doing so, sea level would adjust considerably in response to grounded ice loss due to the small ocean domain; however, these velocities are small (on the order of millimeters per second or less) and do not strongly affect ocean circulation.
In the initial ice state (see Fig. 4(c)), the ice is flowing across the grounding line at 300 m/a, increasing to 700 m/a at the ice shelf front. Initial temperature and salinity are uniformly 1°C and 34.2 psu, respectively. As the flushing time of the cavity is short (on the order of days to weeks) we do not expect this to affect coupled evolution significantly. Movement in the subglacial water layer is far slower so the impact of its initial conditions must be investigated; this is done in Section 6. Fig. 4 summarizes the results of the experiment. Fig. 4(a) shows the evolution of the ice sheet-ice shelf geometry. There is rapid thinning of the ice shelf to a new shape, superimposed on near-steady grounding line retreat. The grounding line retreat occurs due to a thinning of grounded ice, which is a response to an increase of ice flux across the grounding line due to ice shelf thinning and buttressing loss. The movement of a grounding line can be a challenging computational issue in ice sheet modelling (Goldberg et al., 2009), but it has been shown that in the presence of relatively weak beds and narrow channel widths, resolution of 1-2 km is sufficient to represent this process (Gladstone et al., 2012). Thus we are confident that the experiments in this study are able to do so.
The thinning and retreat of the ice sheet is encapsulated in a measure known as Volume above Floatation (VAF), defined as the area integral of  (16). In this example, M min is equal to (900 m) × (1028 kg/m 3 ) and ice thickness decreases linearly with x with a slope of 10 −3 . The − x coordinate is scaled by ε m2 . Δ Mass is mass per unit area relative to M min . (right) The quantity d por (see (20)) as a function of column depth h, with parameter ε por = 1 m. where + (·) indicates the positive part and − (·) the negative part of the quantity. VAF represents the volume of ice that, if melted, would contribute to sea level rise; floating ice makes no contribution. Fig. 4(d) and (f) illustrate this concept: the profiles of actual mass are shown, as well as the effective mass felt by the ocean (cf Fig. 2). The difference between the two profiles is the mass of the ice above floatation. The rate of loss of VAF steadily increases as the ice shelf thins (Fig. 4(b)), but when the grounding line retreats substantially and the ice shelf lengthens, the VAF loss decelerates and eventually slows to zero. While this may be due to increased buttressing from an increasing lateral ice shelf area (e.g., Little et al., 2012), it is more likely because of the diminishing ice sheet having increasingly less volume above floatation to lose.
Melt rate profiles ( Fig. 4(d) and (f)) appear quite different between the early and late stages of the simulation, reflecting the different ice shelf profiles (Fig. 4(c) and (e)). Melt rate is controlled by nearby ocean temperature as well as ocean velocity. Melt rate is therefore higher in deeper waters where water is warmer, and in areas of steeper shelf slope which drive faster flow of melt-buoyed water. This melt-freshened water causes an overturning circulation which is closed by the sponge, but does not strongly affect the stratification in the ocean interior. The jagged appearance of the melt rate is related to the "striped" pattern discussed in Losch (2008), and indicates the transition of the ice shelf base into a new vertical level. The "boundary layer" parameterization of Losch (2008), which averages relevant ocean properties (salinity, temperature, velocity, and surface fluxes) over a full vertical cell thickness, alleviates this variability somewhat in the static shelf case, and Jordan et al. (2017) modifies the parameterization for vertical coupling. This variability would be problematic if it reinforced itself, i.e. if it led to ice shelf thinning in a "sawtooth" pattern which exacerbated the variability. In fact, the opposite is observed to happen. Therefore the stripes do not amplify significantly over time, and importantly do not affect coupled dynamics. It is important to note that despite appearances, the dominant direction of melting is vertically, not horizontally. In reality this is due to the low ice shelf aspect ratio; in MITgcm it is explicit as melt is not applied to vertical cell faces. Grounded ice is lost because the ice shelf, which buttresses the flow, is thinned.
Finally, it is seen that while the under-ice layer of water evolves (Figs. 4(b)), its fluctuations are small, giving confidence that even thinner layers could potentially be used. Here the smallest layer we consider is 1 m.

Variable topography: methodology
The simulation presented in the previous section is a completely synchronous simulation of a marine ice sheet with a dynamic grounding line coupled to an ocean model. It is limiting, however, in that the bathymetry is constant. The true interest in coupled ice sheet-ice shelfocean dynamics arises because of variable topography. In particular, large parts of Antarctica rest on a bed that is well below sea level (more than a km in places) and deepens inland (Fretwell et al., 2013). The effects of buttressing aside, the rate of flow of a marine ice sheet across its grounding line depends strongly on bed depth at that location (Weertman, 1974;Schoof, 2007b). This means that if the grounding line of an ice sheet were to retreat over an inland-deepening bed, the rate of grounded ice loss would increase and so would grounded ice thinning rates, leading to further grounding line retreat. Ice shelf buttressing represents a stabilizing process (Goldberg et al., 2009;Gudmundsson, 2013), but ice-ocean dynamics could potentially affect this stability, leading to increased melt rates as grounding lines retreat and ice shelf cavities change shape (Jacobs et al., 2011;De Rydt et al., 2014).
A difficulty arises, however, in attempting to apply the methodology presented in Section 2 to varying topography. This is because in a zlevel coordinate ocean model, fluid cells can only "communicate" with cells which are directly adjacent. Consider the schematic in Fig. 5. In columns 1 and 2, ice is thicker than its floatation thickness; the subglacial water layer described in Section 2.3 should remain at or near h mwct . In column 3, ice has dropped below its floatation thickness; and in column 4, ice is fully floating. With a flat bathymetry, water would simply flow from the column 4 into the column 3 as the ice thins and the load on the column decreases, and would do so in accordance with the mass and momentum balances of the ocean model. However, there is no hydraulic connectivity between any of the cells in columns 3 and 4, and therefore this cannot happen. As the ice in column 3 thins further, fluid may flow from the column 2, increasing the depth of column 3. There are two possibilities: either column 2 will thin to zero, causing an error; or the fluid cell in column 3 will thicken enough to split into two cells, the upper one of which is hydraulically connected to column 4. However, the pressure differential between columns 3 and 4 could be enormous, leading to tsunami generation and disturbance. Such disturbances are artificial, a result of finite mesh spacing; here we seek to minimize the disturbance.
One solution is to find combinations of horizontal and vertical resolution and subglacial water layer depth that ensure a hydraulic connection between all cells at all times. This could become quite complex, and could lead to either unfeasibly high resolution or subglacial water layers so thick they can no longer be considered negligible. For instance, the bathymetry for a recent ice-ocean coupling intercomparison (Asay-Davis et al., 2016) would lead to layer thicknesses up to 50 m with 1 km horizontal resolution (a high resolution for such a large domain).
An alternative solution, and the one presented in this study, is to allow hydraulic connectivity even between cells that are not at the same z-level. Models of subglacial hydrology (e.g., Schoof, 2010;Werder et al., 2013) follow such a solution, thus avoiding constraints regarding the transfer of fluid between different elevations. They define a hydraulic potential which replaces pressure as a driving force. Our aim is to implement this property within the MITgcm ocean model (in certain cases). We stress that our subglacial water layer is in no sense meant to represent the actual flow of water under ice sheets, which is extremely heterogeneous, with rate factors highly dependent on the direction of flow (Hewitt, 2013). Rather, it is a compromise to enable the use of an ocean GCM for a problem for which it was not originally designed. As such, the implementation of hydraulic fluxes within the layer is simply a means of allowing the layer to fill its purposethat is, to thicken or thin in a locally mass-, salt-, and heat-conservative manner when variations in ice overburden force it to do so.
Returning to the schematic in Fig. 5, in order to allow fluid to move freely from the column 4 to column 3 (as well as columns 1 and 2) when hydraulic potential makes this favorable, several developments are necessary beyond those made in Section 2, and they are discussed below. Fig. 5. Schematic highlighting the difficulty when retreating grounding line over variable topography. In columns 1 and 2 the ice thickness is above the floatation value, whereas in columns 3 and 4 it is below the floatation value. In the columns 1-3 the ocean is masked (there is no fluid) in every level but + k 1, while the rightmost column is masked out in level + k 1. (Despite appearances, there is no fluid in column 3, level k, from a computational perspective; rather, the cell at level k has expanded past its reference surface.) As a result, fluid cannot be transported from column 4 to column 3, even though such transport is favored by the surface load gradient. Fluid can be transported from column 2 to column 3, which is favored by the surface load gradient, but the amount of fluid that can be transported is limited by the depth of column 2.

Porous flux
The issue illustrated by Fig. 5 arises because the ocean floor and ice shelf base are discretized to the ocean grid. In reality, there should in general be pathways and conduits for the water to follow, and if narrow enough these conduits will be controlled by lubricating flow. To represent this flow we modify Eq. (12) to be where Q is a "porous flux", defined at velocity points in MITgcm's Cgrid, equal to where κ por and d por are defined below. Note the second expression above decomposes Φ into its free-surface and baroclinic components, and the third defines ψ low . p′ is the baroclinic component of the pressurethat is, the component due to the vertical integral of density from = z 0. ′ p low is the baroclinic pressure component at the center of the lowest cell of each column. P low is taken from P(z), a horizontal average pressure profile, which in turn is found from a horizontally-averaged density profile ρ z ( ) h (where the h-subscript is to distinguish from the vertical average in Eq 12). ρ z ( ) h is found by averaging over all fluid-filled cells at each depth level (using ρ ref in any levels that contain no such cells). The definition of ψ low may be unfamiliar to readers, and can be reasoned as follows. The potential governing flow in hydrology models (e.g., Schoof, 2010;Werder et al., 2013) can be seen as the deviation of pressure from a reference pressure profileone that is derived via a constant density. Our reference pressure, rather, is derived from a reference density profile, i.e. ρ z ( ) h . Our definition of Φ is used on the basis that its gradients transition smoothly to ocean pressure gradients when the subglacial water layer thickens, as shown below; and that it minimizes movement far inland from the grounding line.
(see also Fig. 2) and κ por and ε por are parameters to be set. This definition ensures the porous flux transitions smoothly to zero as h increases from zero to ε por . The product κ por d por resembles a hydraulic conductivitybut in no way is it meant to represent the true hydraulic conductivity of the subglacial hydrological system. For two adjacent cells at the same z-level, Q is such that if they are "equalised" with respect to the porous flux, there should not be any ocean pressure gradient either. If the grounding line in Fig. 5 retreats further and a connection is established between columns 1 and 2, there will be little to no pressure gradient driving flow. The same is not true for columns 3 and 4, since ψ low corresponds to level + k 1 in column 3, and to level k in column 4. Once connectivity is established (assuming equalisation of the porous flux), the difference in Φ at level k between columns 3 and 4 due to baroclinic effects is where − ρ z ρ z ( ) ( ) h is the deviation from the horizontally-averaged density profile in column 3. With a bedrock step of 50 m over a single grid cell of 1 km width, and an average deviation of ∼ 0.2 kg/m 3 , this will yield a potential difference of ∼ 0.1 m 2 /s 2 , yielding an acceleration of 10 −4 m/s 2 when connection is established. This is not a large acceleration, while it is likely an overestimate of density variations, but this error should be kept in mind.
It is straightforward to implement this porous flux implementation. Eq. (14), the elliptic equation for the updated free surface + η , n 1 becomes The free surface calculation is independent from the thermodynamic component of the ocean model. The above changes enable transport of mass across barriers to the ocean model, but not tracers (temperature and salinity and any passive tracers). We enable this transport with a simple first-order upwind scheme between the two extremal cells on either side of a boundary where porous flux is nonzero, where "extremal" means the exchange is between the top fluidcontaining cell in one column and the bottom of another, or vice versa.
Our parameterization bears similarity to a bottom boundary layer parameterization (BBLp), which has been developed in several variants in ocean models in order to represent small-scale dense overflows (Killworth, 2003). Our parameterization is distinct from most BBLp's in that there is a net transfer of mass as well as tracer. On the other hand, we do not search for "neutral density" levels between which to transfer heat and salt, which is done by some overflow parameterizations (Campin and Goosse, 1999).
The flux Q exists simply because the discretisation of topography prevents a smooth advance of the ocean into new regions opened up by thinning ice. For the rate of grounding line retreat not to be limited by the porous flux parameters κ por or ε por , the flux Q must be large enough that the pressure of the subglacial water layer is roughly equalised with that of the ice shelf cavityi.e., porous flux must not be a limiting factor. It is worth considering typical magnitudes of Q. Assume that, in an ice stream, the upstream region within ∼ 5 km of the grounding line were thinning at an average rate of 50 m/yr (almost an order of magnitude larger than even the fastest-thinning ice streams today). If this region were going afloat, a volume flux of ∼ 0.008 m 2 /s into the newly-created cavity would be required, implying velocities of a similar magnitude if ε por ∼ 1 m. Such velocities would not disrupt the circulation within the cavity. Nevertheless, in subsequent sections we investigate the sensitivity of large-scale behaviour to the porous flux parameters.

Implicit bottom drag
At tracer points, i.e. at cell centers, fluid cell thicknesses are prevented from becoming too thin. However, grounding line migration over variable topography means cell thicknesses at u-and v-points have no lower limit. This is problematic as bottom drag in the MITgcm ocean model is explicit, and this leads to numerical instability in cells that are too thin. A new scheme for implicit bottom drag has been developed in MITgcm. The scheme is a modification of the already-existing scheme for implicit vertical viscosity. As the former is best understood in terms of the latter, an explanation of both is given in the Appendix A.

Variable topography experiment
To illustrate the implementation of the porous flux scheme, an experiment similar to that of Section 3 is carried out, but with variable topography. The bed shallows from −2090 m at = x 0 to −900 m at = x 320 kma slope of 0.005, which is steep for the average bed slope of an ice stream (Fretwell et al., 2013). Ice model input flux is 5 × 10 5 m 2 /yr per unit width. The profile to which the ocean is relaxed is shown in Fig. 3(b), which has slightly cooler water (1°C) at depth, but is otherwise similar. All other parameters are identical to those in Section 3, save those relating to the porous flux parameterization. The experiment is initialized in the same fashion, by first steadying the ice flow model without any basal melting, and then initialising the coupled model with resulting ice thickness.
In Fig. 6, results are shown for a 24-year run with = h 3, mwct κ por = 4.0 s, and ε por = 0.75 m. κ por and ε por are chosen to yield porous fluxes on the order of 10 −4 -10 −3 m 2 /s during fast retreat. Despite the slightly cooler water, the ice stream ungrounds completely in approximately 16 years. The induced circulation is seen to have a much stronger influence on stratification in the cavity than in the flat-bed case. Once the grounding line has retreated and the ice shelf changed shape considerably, multiple overturning cells can be seen (Fig. 6(e)). In this case, the melt plume is reaching neutral buoyancy and flowing out of the cavity at depth, with a distinct plume forming above this depth.
As with the flat-bed experiment, VAF loss rate begins to decrease past a certain point, but again this is likely due to the vanishing of the ice stream. VAF loss rates are much larger than in the flat-bed experiment. Moreover, in the flat bed experiment, the rate of acceleration of VAF loss is constantly decreasing, i.e. the curvature of the plot is positive. This means there is acceleration of VAF loss but it is slowing. Here there is downward curvature around year 7, indicating an instability related to the shape of the bed. As discussed previously, there is known to be an instability related to ice dynamics when a grounding line retreats over a deepening bed. However, it can be seen as well that melt rates increase as the cavity expands, due in turn to an strengthened circulation (Jacobs et al., 2011) and a depression of the freezing point with depth. The subglacial water layer thickness ( Fig. 6(b)) exhibits centimeter-and decimeter-scale fluctuations (rather than the millimeter-scale fluctuations in the flat bed experiment). The short-term fluctuations may be related to ungrounding. If an ice shelf grid cell just downstream of the grounding line is undergoing strong dynamic thinning, it can induce low pressures, which could draw in water from the subglacial water layer. This reduces the layer thickness by no more than a centimeter, however, despite the strong retreat seen in the simulation.
The largest drop in layer thickness is seen when the ice goes completely afloat, which may be due to interaction between the grounding line, the subglacial water layer, and the domain boundary. We anticipate, though, that in normal use the solution will be terminated before the grounding line reaches the edge of the domain. As such, in the sensitivity tests below, we truncate the experiment after 12 years.

Sensitity to porous-flow parameters
A parameterisation has been introduced to an ocean model, and it is not a parameterisation based on a physical process, but rather to circumvent an unrealistic constraint placed upon the ocean by the discretization of topography. It must be shown, then, that the behavior of the model is not sensitive to the parameters involved with the parameterisation; namely κ por and ε por , the hydraulic conductivity and thickness scale of the porous flux layer. The minimum water column thickness h mwct plays a similar role: it is not related to the actual subglacial water thickness but exists instead to prevent the water column from becoming negative. We test the impact of all three parameters on the outcome of the experiment of Section 5. As discussed in Section 4.1, the influence of κ por and ε por should be minimal as long as fluxes are (a) large enough to keep pace with the thinning of grounded ice but (b) small enough that any induced movement is low. Table 2 details the values investigated in experiments 1A-F. Each parameter is varied independently: κ por is halved, and ε por is doubled, relative to Section 5, and h mwct is decreased to 2 m and to 1 m. This parameter variation is seen to affect VAF by at most 0.15% over 12 years of simulation ( Fig. 7(a)). Melt rates vary by at most ∼ 1% (Fig. 7(b)). At the end of the 12 years, during which the grounding line has retreated ∼ 80 km and ice has thinned over a kilometer in places, the variation leads to differences of at most 6 m in ice thickness (Fig. 7(c)). Overall, these purely numerical parameters seem to have very little effect on the evolution of the coupled system. The reference density ρ ref is also varied. A small change in this parameter has little effect on ocean circulation, but as discussed above it affects hydraulic potential gradients in the subglacial water layer. The effect of decreasing ρ ref by 3 kg/m 3 is quite small. Notably, it thins the subglacial water layer slightly close to the upstream boundary (Fig. 7(c)).
It is important to establish that coupled evolution does not depend critically on melt rates in the near vicinity of the grounding line, as this is where circulation is most likely to be affected by the presence of the subglacial water layer and porous flux parameterisation. To this end we mask out melting when the ocean column is below a certain thickness (z melt ), and vary this threshold between 0 and 10 m (see experiments 2A-C in Table 2). The effect on coupled evolution is very small (Fig. 7(d)-(f)). A number of studies have suggested that, in the absence of tidal mixing or significant subglacial runoff, melt rates are small within a short distance (3-5 km) of the grounding line, as water flowing up the underside of the shelf must gain sufficient buoyancy before substantial mixing of heat to the ice-ocean interface can occur Seroussi et al., 2017). Thus we are confident that this relative unimportance of melt rates where the ocean column is < 10 m will carry to other settings. Finally, the initial conditions of the subglacial water layer are investigated; as mentioned previously, movement in the layer is negligible far from the grounding line, as is transport of heat and salt; therefore initial conditions in the layer will persist. We carry out two additional sensitivity tests (experiments 3A-B in Table 2). In these experiments initial temperature in the layer either linearly increases to 1.0°C, or decreases linearly to −0.8°C, with distance from the grounding line. The range of variability seen is similar to that of other sensitivity experiments (Fig. 7(d)-(f)). A possible reason is that where the ocean column is thin (either under grounded ice or very close to it), movement is so slow that thermodynamic processes, such as diffusion mixing and very low rates of melt, are sufficient to homogenize the ocean conditions and shield the cavity from the conditions of the subglacial water layer. The impact of not using the porous flux parameterization at all can be observed by setting all relevant parameters to zero (except h mwct , which is necessarily nonzero to allow grounding line retreat, see experiment 4 in Table 2). The results are shown in Fig. 8. The result is shown from the perspective of the ocean, for which the grounding line does not retreat at all despite significant thinning of the ice shelf. This can be understood in reference to Fig. 5: due to ice shelf thinning, the ice in column 3 has thinned sufficiently to go afloat; but for the ocean column to inflate, water must flow in from column 4, and it cannot. From the ice model perspective, the grounding  line has retreated. As mentioned in Section 2.2, the ice model calculates its own grounding line irrespective of the ocean model; the ocean model only influences the ice model through melting. (As plotted in Fig. 8(a), the ice surface appears to have a local maximum, but this is only from the ocean perspective.) Nevertheless, the ice sheet has not thinned as much as in Fig. 6(e), as the ocean cavity has not increased in extent. If the ice model were to determine its grounding line position based on ocean column thickness rather than ice thickness, there would still be little grounding line retreat and the ice would have thinned even less; but as mentioned in Section 2.2, this modification is not made to the ice model in the present study. We infer that, at least for our simplified two-dimensional experiment, coupled evolution is insensitive to porous flux parameters within a reasonable range; but excluding porous flux altogether over steep topography can be problematic.

Advancing grounding line
In all prior experiments, coupled grounding line retreat has been demonstrated; but for our method to be robust, it must be shown that grounding line advance can be simulated as well. We carry out three additional experiments where the steep-topography runs detailed above are halted after 8 years and exposed to a step-change in far-field oceanic conditions by lowering temperature to −1.9°C uniformly at the boundary.
One additional change is made: the buttressing half-width W is decreased from 25 to 15 km. The reason for this change is that it was found the ice stream has reached a point of instability after 8 years of retreat. From this point, the grounding line will continue its retreat without stopping, even if melt rates are set to zero. This instability is often referred to as the Marine Ice Sheet Instability (Weertman, 1974;Schoof, 2007a), in which grounding line retreat on a seaward-sloping bed accelerates as the grounding line deepens. As much of the Antarctic margin has similarly-sloped beds, the instability is cause for concern. However, our aim is not to verify this instability but capture coupled evolution in the presence of a retreating and advancing grounding line; thus the step-change in W, while artificial, is necessary to cause grounding line advance.
In contrast to the high rate of retreat in the previous experiments, advance occurs more slowly, advancing ∼ 10 km in 16 years ( Fig. 9(a)), but is sufficient to establish that advance can be represented. A full parameter sensitivity test is not carried out but we examine sensitivity to h mwct (Fig. 9(c)-(e)). As the primary driver of change in VAF is the step-change in W, the sensitivity of VAF to layer thickness might not be indicative of the robustness of the algorithm. Melt rate and ice shelf draft, however, are insensitive to layer thickness, just as in the retreat experiments.

Three dimensional coupled experiment
A three-dimensional coupled experiment is carried out over variable topography. The domain is 300 km long by 50 km wide with 1 km horizontal resolution in ice and ocean models, and 20 m vertical resolution in the ocean model. Bathymetry is similar to that described in Section 5. It is uniform in the across-flow (x) direction and rises in the ydirection with a slope of 0.005 from -2000 m to −900 m over 220 km, the location of the ice front. It then levels off, and beyond this point open ocean is assumed. Ice enters the domain at the southern boundary at a uniform flux of 2.25 × 10 6 m 2 /yr per unit width. The profiles of temperature and salinity to which the ocean is relaxed in the sponge layer are as in Fig. 3, with the same restoring time constant. The ice model has no-slip conditions along the x-boundaries. In the ocean, the minimum water column thickness (h mwct ) is 1.0 m, the porous flux conductivity κ por is 4 s, and the porous flux layer thickness ε por is 0.75 m. The experiment is carried out on the f-plane with Coriolis parameter f set to − − 10 4 . All other ice and ocean parameters are as described previously or in Table 1.
The results are presented in Fig. 10. Overall evolution is similar to the two-dimensional experiment, though initially a "bulge" forms at the base of the ice shelf and propagates to the front ( Fig. 10(a)), having a small but noticeable effect on local melt rates due to its effects on subice velocity. It is not known whether such features form in real ice shelves as they are transient in our model, though similar features were observed in the experiments of Jordan et al. (2017), despite there being no moving grounding line. After 40 years, the ice stream is nearly   Goldberg et al. Ocean Modelling 125 (2018) 45-60 completely afloat. Melt rate response is similar to the two-dimensional runs: after an initial ∼ 5-year adjustment the bulk melt rate increases ( Fig. 10(b)), which coincides with the onset of grounding line retreat and cavity widening/deepening. VAF loss rates are considerable -50 km 3 /a, which would translate to ∼ 0.14 mm/a global sea level rise though this number should not be seen as representative, as the ice shelf was initially in equilibrium with no melting at all. Grounded ice loss rates lessen toward the end of the runbut as with previous experiments, this is likely due to the grounding line approaching the upstream boundary. There is a strong overturning circulation under the ice shelf which strongly alters the stratification (Fig. 10(c) and (e)), similarly to the two-dimensional experimentsthough this occurs primarily at the western boundary (x = 0) due to the effects of rotation. The Coriolis effect can most easily be seen in under-ice flow and melt rate ( Fig. 10(d) and (f)). Most melting takes place at depth where water is warm, but as geostrophy prevents the buoyant melt-laden water from flowing directly up the sloped ice shelf base it flows to the western edge of the domain, where boundary effects break the geostrophic constraints. As a result, velocities are quite high in this boundary, as are melt rates, yielding considerably more thinning at the western boundary than elsewhere in the shelf. Ocean-driven melting has a strong effect on ice shelf flow as well, as can be seen in Fig. 11. The initial ice configuration (i.e. in steady-state with the absence of melting) yields a symmetric velocity which moves fastest at the ice shelf front and has very little flow in the transverse direction. After 40 years, the flow pattern is strongly asymmetric, with speeds increasing toward the western side, and flow is significantly deflected in the western direction. These features can be largely explained by the melt-driven channel along the western boundary, which allows greater shear within the ice shelf and also induces a flow pattern which acts against the thickness differential. This asymmetry, however, does not express itself in the flow pattern across the grounding line, which explains why the grounding line remains roughly straight throughout the simulation.

Discussion
The major output of this study is a coupled ice sheet-ocean model that can represent grounding line migration over variable topography, for which the ocean component can run continuously without requiring mass, salt or heat to be arbitrarily moved to different locations. Numerical experiments were carried out in which far-field forcing was applied to the ocean model, inducing retreat of an ice sheet over a bathymetry which deepens inland. It is fair to point out that these simulations could have been carried out with an asynchronous modeland, in fact, similar experiments have been carried out, yielding similar qualitative behaviour. However, the locally conservative properties of our scheme additionally mean that our methodology can be implemented in a regional or global ocean model without requiring artificial sources and sinks to preserve conservationwhich would be required with an asynchronous approach. Additionally, asynchronous approaches still leave open the question as to whether ocean dynamics are affected by the intermittent coupling between ice and ocean. In effect, our approach reduces the coupling time step to the ocean time stepand all approaches with longer coupling time steps are an approximation to this solution. Moreover, in the experiments presented here, ocean forcing was persistent. But in reality forcing of various time scales does occur, from centennial time scales reflecting climate changes (Boning et al., 2008) to monthly or potentially weekly time scales reflecting fast glacial processes and subseasonal ocean time scales (Smith et al., 2017;Christianson et al., 2016). With an asynchronous model, forcing the ocean model on time scales on the order of, or shorter than, the coupling frequency would be problematic as aspects of the coupled evolution would not be captured. Our synchronous coupling allows for forcing on daily or even hourly time scales. Our mode of coupling is not yet complete, however. Marine ice sheet models determine the location of the grounding line, and the depth of the ice shelf draft, by assuming a static ocean of uniform density; the same is done in our ice sheet model, which is influenced by the ocean through melting only. A future plan for development is to allow the floatation of ice to be determined by the ocean state.
A further advantage of a synchronously coupled ice-ocean model is the prospect of an adjoint coupled model. Adjoint models allow detailed investigation of sensitivitiesin a manner that would be intractable through standard parameter variationand enable state estimation and uncertainty quantification. Adjoints of both ice and ocean components are easily and flexibly generated through algorithmic differentiation (Heimbach and Losch, 2012;Goldberg and Heimbach, 2013). Synchronous coupling between the two now allows for a coupled adjoint. A coupled adjoint would allow detailed sensitivities of the coupled iceocean system to be investigated, and could aid in initialization of the coupled system to observations. In our implementation a thin layer of ocean is maintained under grounded ice. This is essentially a numerical approximation so that the ocean can advance continuously into previously ice-covered groundand as with all approximations, the thickness of the layer should not strongly impact the model evolution. In our sensitivity study in Section 6 we vary the thickness of this layer between 1 and 3 m and see little impact. In order to allow movement within the layer in the Fig. 10. Output of the three-dimensional coupled experiment on sloped bathymetry. (a) Horizontally-averaged ice sheet geometry at annual intervals.
Color indicates year of experiment. (b) Time series of integrated melt rate and VAF, as well as the minimum thickness of the ocean column. (c) Ocean salinity (shading) and overturning streamfunction (black contours; spacing 10 5 m 3 /s) at the western boundary, after 10 years of simulation. Maximum velocity is 1.28 m/s. (d) Melt rate (shading), under-ice velocity (vectors), and ice shelf draft (black contours). Geometry is presented as deviation from the width-averaged draft: thick contour is 0, with thin solid contours spaced 20 m (up to 100m) for positive deviation and dashed contours spaced 20 m (down to −100 m) for negative deviation. (e,f): same as (c,d) at 40 years. Maximum ocean velocity is 1.2 m/s. The "sawtooth" pattern of the overturning stream function is due to partial cells not being taken into account in the diagnostic calculation, but this has no effect on the evolution of the model. presence of variable topography, we introduce a porous flux parameterisation that is driven by a potential which is related to ocean pressure. Without this parameterisation, the layer would need to be considerably thicker than the values stated above in the presence of steep topography. The active ocean layer under grounded ice should not be confused with a true subglacial water layer, which is fed by basal melting of the ice sheet. Still, the movement of such water is driven by a balance between friction, overburden, and potential gradientwhich is essentially what drives the water in our thin ocean layer. As such it may be possible to modify the dynamics of the layer slightly to properly represent a subglacial hydrological regime which continuously transitions to an ocean cavity regime. This would require significant further investigation and development effort, however.

Conclusions
We develop a scheme that allows for synchronous coupling between a marine ice sheet model with a dynamic grounding line and a z-coordinate ocean model. While "vertical coupling" (the continuous evolution of the ice shelf/ocean interface due to melting) was developed in a separate study , this study focuses on "horizontal coupling", i.e. the advance and retreat of the ocean's lateral extent.
With a flat bathymetry, little modification to the ocean model is required (beyond limiting ice overburden to allow a thin subglacial water layer). In the presence of variable topography, a "porous flux" parameterization is required to allow the ocean to infiltrate new areas where ice has thinned to the point of floatation. The porous flux is not based on a physical model; rather, it exists to counter the nonphysical obstruction of water due to the discretization of topography. The porous flux transport is small enough that it has negligible effect on the ocean circulation. Yet, it is large enough that ice dynamics remains the controlling factor of grounding line migration in realistic simulations.
The modifications made to the code are relatively non-invasive: the major modification is to the Poisson equation solved for the implicit free surface evolution. This introduces a few unconstrained parameters, but sensitivity of model behaviour to these parameters is shown to be small or negligible. Experiments carried out in two and three dimensions demonstrate physically reasonable behavior, suggesting the coupled model is suitable for scientific experimentation.
In the ocean model L is discretized using a finite volume approach, yielding a tridiagonal matrix for each column that can be easily inverted independently (e.g. using L.
where the second equality is due to the particular form of L which only acts on vertical shear without changing the depth integrated transport.

A.2. With implicit surface or bottom drag
Here we describe our modification to the above scheme to allow for implicit bottom drag. (23) is rewritten as D.N. Goldberg et al. Ocean Modelling 125 (2018) L d can be discretized similarly to L, the one difference being the addition of the surface and/or bottom drag terms to the diagonal matrix coefficient in the relevant rows. The inversion of L d in the second term on the right hand side of (28) no longer reduces to the identity but instead (28) can be written as where I is a function equal to unity in the vertical. Note L − I { } d 1 is a function of depth which we will refer to as J.
where J n is the vertical average of J at time n, which becomes a weighting factor in the elliptic operator that is solved for the free-surface update.