BISMG:VickyL

Flow beneath ice shelves
Aim is to model the flow of ocean water into the cavity beneath an ice shelf and its interaction with the overlying ice. We want to capture the effect of intrusions of warm water from the surrounding oceans on the basal melt rate of the ice shelf.

Conceptual model
Melt water from the base of the ice shelf and from the foot of an ice stream is assumed to form a layer adjacent the underside of the shelf. The layer flows along the underside to the shelf's front. As it flows ocean water in the cavity is entrained into the melt-water layer.

Model Assumptions

 * 1) Pressure is assumed hydrostatic.
 * 2) Pressure is continuous across layer interfaces

Mathematical model

 * 1) $$\nabla . \Big( \frac{1}{D+H} \nabla \psi \Big) = \zeta$$
 * 2) $$\zeta = - \frac{f}{c_{d}} \nabla . \overline{\boldsymbol{u}} - \frac{1}{\rho_0 c_{d}} \boldsymbol{k} . \nabla \times \overline{\nabla p_{b}}$$
 * 3) $$\overline{u} = - \frac{1}{(D+H)}\frac{\partial \psi}{\partial y} \qquad \overline{v} =  \frac{1}{(D+H)}\frac{\partial \psi}{\partial x}$$
 * 4) $$f \boldsymbol{k} \times \bar{\boldsymbol{u}_{p}} + c_{d} \bar{\boldsymbol{u}_{p}} =   f \boldsymbol{k} \times \overline{\boldsymbol{u}}  +  c_{d} \overline{\boldsymbol{u}} +   \frac{g}{ \rho_{0} (D+H)} \left[ \frac{DH}{2} \nabla \rho_{p} - \right.  \left. H (\rho_{c} - \rho_{p}) \nabla I + \frac{H^{2}}{2} \nabla \rho_{c} \right] $$
 * 5) $$\bar{\boldsymbol{u}}_{c} = \frac{D+H}{H} \overline{\boldsymbol{u}} - \frac{D}{H} \bar{\boldsymbol{u}}_{p}$$
 * 6) $$\frac{\partial D}{\partial t} + \nabla . \left( D \bar{\boldsymbol{u}_{p}} \right) = \dot{e}$$
 * 7) $$H = h - B - D$$
 * 8) $$\frac{\partial T_{p}}{\partial t} + \bar{\boldsymbol{u}_{p}} \nabla . T_{p} = K_{h} \nabla^{2} T_{p} + \frac{\gamma_{T}|\bar{\boldsymbol{u}_{p}}|}{D} (T_{B}-T_{p})  + \frac{\dot{e}}{D} (T_{I}-T_{p})$$
 * 9) $$\frac{\partial S_{p}}{\partial t} + \bar{\boldsymbol{u}_{p}} \nabla . S_{p} = K_{h} \nabla^{2} S_{p} + \frac{\gamma_{S}|\bar{\boldsymbol{u}_{p}}|}{D} (S_{B}-S_{p})  + \frac{\dot{e}}{D} (S_{I}-S_{p})$$
 * 10) $$\frac{\partial T_{c}}{\partial t} + \bar{\boldsymbol{u}_{c}} \nabla . T_{c} = K_{h} \nabla^{2} T_{c} + \delta_{D} \frac{\gamma_{T}|\bar{\boldsymbol{u}_{c}}|}{H} (T_{B}-T_{c})  - \frac{\dot{e}}{H} (T_{I}-T_{c})$$
 * 11) $$\frac{\partial S_{c}}{\partial t} + \bar{\boldsymbol{u}_{c}} \nabla . S_{c} = K_{h} \nabla^{2} S_{c} + \delta_{D} \frac{\gamma_{S}|\bar{\boldsymbol{u}_{c}}|}{H} (S_{B}-S_{c})  - \frac{\dot{e}}{H} (S_{I}-S_{c})$$

Indexing
PetSc routines are written in C where vectors and arrays have a lower index of 0. Fortran 90 vectors and arrays start from 1 by default. This how to put a Petsc Vec variable into a 1d allocatable array

real(dp), dimension, allocatable :: f90arr Vec :: petscvector PetscInt :: N PetscInt, dimension, allocatable :: indx

! Create the PetSc vector call VecCreateSeq(PETSC_COMM_SELF,N,petscvector,ierr)

... fill petscvector

allocate(indx(N), f90arr(N))

! Create array of indices for petscvector to extract the data do i=1,N k=i-1 indx(i)=k end do

! Put data in petscvector into the f90 array f90arr call VecGetValues(petscvector,N,indx,f90arr,ierr)

f90arr(1) contains the first element of petscvector f90arr(0) is undefined

Work in progress
Suggestions and comments from JCRP meeting (13/11/09):
 * Why is diffusion of scalars (T,S) important near the inflow along the grounding line? Surely the scalar difference is driving the plume?
 * Use the simplified momentum equation derived under the assumption of stationary cavity for the motion of the plume with the moving cavity. We know this equation works. Continuity of the whole cavity is satisfied through the streamfunction. Velocity in the cavity is simply the difference between the depth weighted velocity in the whole cavity calculated from the streamfunction and the plume velocity. Question: can we really justify using a stationary ocean assumption just to calculate the plume dynamics and then ignore it?
 * Why not ignore the pressure gradient due to the base of the ice shelf and force term due to the slope of the ice shelf in plume momentum? These terms exist even if there is no flow beneath the ice shelf - they do not drive motion so they should be ignored. In the above model the pressure gradient and the force due to the slope of the ice shelf is found from the barotropic momentum equation and then substituted into the plume equation. Crucially, this substitution introduces a coupling to the cavity flow through the reduced gravity term. By ignoring the two 'forces' the cavity only influences the plume through entrainment in the thickness and scalar equations. There is a reduced gravity term in the momentum equation for the cavity. Perhaps we could calculate the cavity flow directly and then the plume velocity would be the difference between the depth weighted barotropic and cavity velocities.
 * What would happen if the cavity flow went to zero? I'm not sure how to switch off the cavity flow in the above model. It was suggested that the drag coefficient for the cavity be increased to slow the flow down. The cavity momentum is not calculated directly in this model. I would need to increase the drag in the barotopic momentum equation.

The week ahead

 * Need to distinguish between land/sea boundary and open boundary when differentiating. Do this through masks. How about inflow boundary on the grounding line?
 * Investigate the inflow boundary on the grounding line.
 * Play with solvers in PetSc subroutine.
 * Implement extrapolated open boundary conditions.
 * Try solving the plume momentum equations in the form of linear equations by introducing $$\frac{ \partial \boldsymbol{u_{p}}}{ \partial t}$$ and marching the solution forward rather than treating the momentum equations as two simultaneous equations. This has the advantage of not requiring cross derivatives.
 * Spurious oscillation:
 * 1) try upwinding rather than central differing
 * 2) try radiative open boundary conditions. A mixture of the field variable at different times.
 * 3) try solving the momentum equations as above and introduce nonlinear terms.

Friday 11th December 2009

 * Upwinding didn't stop oscillation. Delayed onset as the velocity and depth were damped slightly.
 * Extrapolation didn't work for PIG. Got negative depths and spurious velocities. Wave developed across the boundary.
 * North and east boundaries handled slightly differently from south and west ones.
 * Tried radiative boundary condition. Seems to generate an inflow and plume depth increases at boundary.

Friday 4th December 2009
$$f \boldsymbol{k} \times \bar{\boldsymbol{u}_{p}} + c_{d} \bar{\boldsymbol{u}_{p}} =  f \boldsymbol{k} \times \bar{\boldsymbol{u}_{c}}  +  c_{d} \bar{\boldsymbol{u}_{c}} +   \frac{g}{ \rho_{0}} \left[ \frac{D}{2} \nabla \rho_{p} - \right. \left. (\rho_{c} - \rho_{p}) \nabla I + \frac{H}{2} \nabla \rho_{c} \right] $$
 * Modified derivatives of scalars onto scalar grid points, i.e. gradients used in Jacobians, to take account of known boundary conditions. T, S and $$\rho$$ are insulated along the land/sea boundary so the normal gradient is zero. The normal derivative of the scalar adjacent to wall is calculated by first calculating the gradient across the cell face opposite the wall using central difference and then taking the mean value (or linear interpolating) with the known value on the boundary.
 * The momentum equation for plume flow above a stationary, homogeneous cavity is recovered by substituting an expression for barotropic velocity in terms of plume and cavity velocities (eqn 5) into the momentum equation (4) above

and setting $$\bar{\boldsymbol{u}_{c}}=0$$ and $$\nabla \rho_{c}=0$$
 * Tested the plume velocity by switching off the cavity coupling as above on the simple geometry set up (bathymetry type 12). Switching off the couple made very little difference to the flow pattern - the depth of the plume was still shallow and did not travel along the western boundary and the velocity field was weak compared to the reduced physics model. The difference between the current and reduced physics model is the drag term. The current model drag uses a linear law $$c_{d} \bar{\boldsymbol{u}_{p}}$$ but the reduced physics model uses a quadratic law and is a function of depth $$ \frac{c_{d} |\bar{\boldsymbol{u}_{p}}|}{D} \bar{\boldsymbol{u}_{p}}$$. For the reduced model, plume speed is about ~0.1 m/s or less and depth is ~100 m along the western boundary and ~10 m elsewhere, so the magnitude of the drag term here is much smaller than that of the current model. This accounts for the difference between the two flow patterns.
 * Reduced the drag coefficient for the current model with the simple geometry set up. Reducing it by a factor 10 to 0.0003 was okay, though a little messy around the inflow. Much smaller and the code crashed. Depth and velocity of the plume did not greatly increase.
 * Went back to the reduced physics model and looked at solving the velocity fields separately rather than using the solution to the simultaneous equations. One can arrange the equations such that either the velocity in the drag term is the subject or in the Coriolis force. Using the drag term means that only the along stream derivatives are required. However, since the drag is quadratic we would have to divide by the linearised speed which rules this approach out (initial velocity is zero). Using the Coriolis force only the cross derivatives are used. It was particularly problematic near the left hand corner of the inflow for the simplified geometry set up (bathymetry type 12). Numerical oscillation appeared to radiate away from the wall. When the solution to the simultaneous equation is used this is not a problem. The use of derivatives in two directions seems to (spatially) stabilise the flow. It should be noted that even with this method the velocity at the inflow alternates values with each timestep (i.e. temporal effect) for the initial 'spin up' and then settles down.
 * Noticed problems with the reduced physics model for the simplified geometry set up (bathymetry type 12) at long time. The open boundary conditions are not very good; the flow is not leaving the computational domain probably. Instead the boundary is behaving like a porous wall. Also, a spurious oscillation appears to start on the top left hand corner where a wall meets the open boundary and travels along the left hand wall. The plume is relatively fast moving and deep here. This oscillation is large but does not destabilise the flow. It is not present in the fully nonlinear model.
 * Replaced the no normal gradient open boundary condition for the northern boundary of the simplified geometry set up (bathymetry type 12) with extrapolation of T, S, $$\boldsymbol{u}_{p}$$ and D. Much improved. Flow is leaving domain properly. Need to alter the implicit solver to include extrapolation.
 * Spurious oscillation still a problem.

Friday 27th November 2009

 * Played with a PetSc example of a general elliptical problem, the 2d diffusion equation. Wrote a code to solve the example on p261 of Numerical Solution of PDEs, Morton & Mayers. Main difficulty was extracting the solution from a PetSc variable to a fortran array. Stephen helped.
 * Included a petsc subroutine into the plume model to solve the streamfunction-vorticity equation. Initial runs with PIG and simple geometry (bathymetry_type = 12) set ups suggest petsc works with Dirichlet and Neumann boundary conditions.

Friday 20th November 2009

 * Slapsolv cannot solve the streamfunction-vorticity equation with a Dirichlet condition on the land/sea boundary AND a Neumann condition on the open boundary. Peric's linear solves can handle the problem but are painfully slow. Need to try something else.
 * Gethin installed PetSc on Dartagnan.
 * Read through PetSc manual and selected a suitable example.

Friday 13th November 2009

 * Adjusted the cross derivatives of scalar quantities onto velocity grid points. Previously, I used second order one-sided finite difference using 3 scalar points adjacent to the boundary or first order one-sided difference if only 2 scalar point were available. Now the derivatives use known boundary values if available. For the streamfunction the value on the land boundary is known so I use second order difference formulation for a nonuniform grid (the spacing between the boundary value and the adjacent cell centre is $$\Delta x/2$$. For T,S and $$\rho$$ the normal gradient to the land is zero (i.e. insulating boundary). I calculate the gradient half a grid cell above the required velocity point using central difference and then average it with the boundary value (i.e. divide by 2 since the boundary value is zero). The results were not that different from the previous version. It did reduce the variations the divergence of barotropic velocity used in the streamfunction-vorticity equation around the grounding line in the PIG simulations.
 * Introduced two-way entrainment.

$$ \dot{e} = \begin{cases} \frac{c_{L}^{2}}{S_{m}} \left( |\bar{\boldsymbol{u}_{p}}-\bar{\boldsymbol{u}_{c} }|^{2} + \frac{g ( \rho_{c}-\rho_{p} ) D}{\rho_{0} S_{m}}\right)^{\frac{1}{2}} & if \,\, D > D_{\min} \,\, and \,\, |\bar{\boldsymbol{u}_{p}}| > |\bar{\boldsymbol {u}_{c}}| \\ -\frac{c_{L}^{2}}{S_{m}} \left( |\bar{\boldsymbol{u}_{p}}-\bar{\boldsymbol{u}_{c} }|^{2} + \frac{g ( \rho_{c}-\rho_{p} ) H}{\rho_{0} S_{m}}\right)^{\frac{1}{2}} & if \,\, D > D_{\min} \,\, and \,\, |\bar{\boldsymbol{u}_{p}}| \le |\bar{\boldsymbol {u}_{c}}| \\ 0 & otherwise \end{cases} $$

where the Schmidt number is a function of the Richardson number which is now defined as

$$ Ri = \begin{cases} \max \left( 0.05, \frac{g ( \rho_{c}-\rho_{p} ) D}{\rho_{0} |\bar{\boldsymbol{u}_{p}}-\boldsymbol{u}_{c}|^{2}} \right) & if \,\, |\bar{\boldsymbol{u}_{p}}| > |\bar{\boldsymbol {u}_{c}}| \\ \max \left( 0.05, \frac{g ( \rho_{c}-\rho_{p} ) H}{\rho_{0} |\bar{\boldsymbol{u}_{p}}-\boldsymbol{u}_{c}|^{2}} \right) & if \,\, |\bar{\boldsymbol{u}_{p}}| \le |\bar{\boldsymbol{u}_{c}}| \end{cases} $$
 * Values of T and S at the interface are set to the depth weighted average of the plume and cavity values such as $$T_{I}= \frac{D T_{p} + H T_{c}}{D+H}$$. Previously use the mean value between the plume and cavity values. Result for PIG is that temperature in the cavity near the inflow on the grounding line does not drift above its maximum initial value (i.e. 1 degree C).