Difference between revisions of "Blatter-Pattyn model"

From Interactive System for Ice sheet Simulation
Jump to: navigation, search
 
(34 intermediate revisions by 3 users not shown)
Line 1: Line 1:
The starting point for the Blatter-Pattyn model is the full Stokes equations
+
The higher-order dynamics scheme implemented in Glimmer/CISM is based on the so-called "Blatter-Pattyn" model. The starting point for this model is the full Stokes equations
  
  
Line 18: Line 18:
  
  
There are a number of ways to argue that because of the "shallowness" of ice sheets - that is because the ratio of <big>''H''/''L''</big>, where ''H'' is the thickness and ''L'' is a relevant horizontal length scale, is small - the equations above can be reduced to the following "first-order approximation"
+
There are a number of ways to argue that, due to the "shallowness" of ice sheets (because the ratio of ''H''/''L'', where ''H'' is the thickness and ''L'' is a relevant horizontal length scale, is small), the equations above can be reduced to the following "first-order approximation":
 
+
  
  
Line 29: Line 28:
  
  
 
+
The arguments that support this reduction are fairly complex, either based on a variational analysis or an asymptotic analysis. Two of the papers given in the references below ([http://www.ingentaconnect.com/content/igsoc/jog/2010/00000056/00000197/art00010 ''Dukowicz et al.'',2010] and [http://qjmam.oxfordjournals.org/content/63/1/73.abstract ''Schoof and Hindmarsh'',2010], respectively) provide more details on the mathematical background that allows us to state that these equations are "first order accurate" approximations to the Stokes equations.  
Unfortunately, the arguments that support this reduction are complex and we can't go into them in detail here. Two papers given in the references below (''Schoof and Hindmarsh'' (in press) and ''Dukowicz et al.'' (submitted)) provide more details on the mathematical background that allows us to state that these equations are "first order accurate" approximations to the Stokes equations.  
+
  
  
Line 39: Line 37:
  
  
This expression can be substituted into the horizontal pressure gradient terms above to remove pressure from the equations. For example, for the 'x' component of velocity we have
+
This is simply a statement the the full vertical normal stress is balanced by the hydrostatic pressure (the so-called ''hydrostatic assumption''). This expression can be substituted into the horizontal pressure gradient terms above to remove pressure from the equations. For example, for the ''x'' component of velocity we have
  
  
Line 51: Line 49:
  
 
Now we can use the incompressibility and stress-strain colinearity constraints to rewrite the vertical normal deviatoric stress in terms of horizontal normal deviatoric stresses. This allows us to write the horizontal balance equations in terms of horizontal stress gradients only. Again, for the ''x'' direction we have
 
Now we can use the incompressibility and stress-strain colinearity constraints to rewrite the vertical normal deviatoric stress in terms of horizontal normal deviatoric stresses. This allows us to write the horizontal balance equations in terms of horizontal stress gradients only. Again, for the ''x'' direction we have
 +
 +
 +
<math>\begin{align}
 +
  & \tau _{zz}=-\tau _{xx}-\tau _{yy} \\
 +
& -\frac{\partial \tau _{zz}}{\partial x}=-\frac{\partial }{\partial x}\left( -\tau _{xx}-\tau _{yy} \right)=\frac{\partial \tau _{xx}}{\partial x}+\frac{\partial \tau _{yy}}{\partial x} \\
 +
& \frac{\partial \tau _{xx}}{\partial x}-\frac{\partial \tau _{zz}}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=\rho g\frac{\partial s}{\partial x}\quad ...\text{becomes}... \\
 +
& 2\frac{\partial \tau _{xx}}{\partial x}+\frac{\partial \tau _{yy}}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=\rho g\frac{\partial s}{\partial x}\quad  \\
 +
\end{align}</math>
 +
 +
 +
Note that at this point we have removed the vertical balance equation entirely; it has been incorporated into the horizontal balance equations using incompressibility. The horizontal balance equations are now:
 +
 +
 +
<math>\begin{align}
 +
  & x:\quad 2\frac{\partial \tau _{xx}}{\partial x}+\frac{\partial \tau _{yy}}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=\rho g\frac{\partial s}{\partial x}\quad  \\
 +
& y:\quad 2\frac{\partial \tau _{yy}}{\partial y}+\frac{\partial \tau _{xx}}{\partial y}+\frac{\partial \tau _{xy}}{\partial x}+\frac{\partial \tau _{yz}}{\partial z}=\rho g\frac{\partial s}{\partial y} \\
 +
\end{align}</math>
 +
 +
 +
Next, we want to write these equations in terms of velocities, since that is what we are ultimately solving for. The link between stress gradients and velocities is through the constitutive law for ice (here we'll assume Glen's law) - relating strain rates to stresses - and the definition of the strain rate tensor - relating strain rates to velocity gradients.
 +
 +
 +
 +
<math>\begin{align}
 +
  & 1.\quad \tau _{ij}=B\dot{\varepsilon }_{e}^{\frac{1-n}{n}}\dot{\varepsilon }_{ij},\quad B=B(T) \\
 +
& 2.\quad \dot{\varepsilon }_{ij}=\frac{1}{2}\left( \frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}} \right) \\
 +
& 3.\quad 2\dot{\varepsilon }_{e}=\dot{\varepsilon }_{ij}\dot{\varepsilon }_{ij} \\
 +
& 4.\quad \eta \equiv \frac{1}{2}B\dot{\varepsilon }_{e}^{\frac{1-n}{n}} \\
 +
& 5.\quad \tau _{ij}=2\eta \dot{\varepsilon }_{ij} \\
 +
\end{align}</math>
 +
 +
 +
 +
In order, the four expressions above give
 +
 +
# Glen's flow law (actually, the inverse flow-law here) where ''B(T)'' is a temperature dependent pre-factor
 +
# The definition of the strain rate tensor in terms of velocity gradients
 +
# The definition of the effective strain rate, <math>\dot{\varepsilon }_{e}</math>, a norm of the strain rate tensor
 +
# A definition for the "effective viscosity" (after rearranging some terms in (1))
 +
# ... which allows us to write the relationship between stress and strain in a standard "Newtonian" way (but with a non-Newtonian viscosity)
 +
 +
 +
Now taking these definitions into our stress balance equations above and expanding in terms of strain rates, velocity gradients, and effective viscosities, we have (for the ''x'' direction only)
 +
 +
 +
<math>\begin{align}
 +
  & x:\quad 2\frac{\partial \tau _{xx}}{\partial x}+\frac{\partial \tau _{yy}}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=\rho g\frac{\partial s}{\partial x} \\
 +
& x:\quad 2\frac{\partial }{\partial x}\left( 2\eta \dot{\varepsilon }_{xx} \right)+\frac{\partial }{\partial x}\left( 2\eta \dot{\varepsilon }_{yy} \right)+\frac{\partial }{\partial y}\left( 2\eta \dot{\varepsilon }_{xy} \right)+\frac{\partial }{\partial z}\left( 2\eta \dot{\varepsilon }_{xz} \right)=\rho g\frac{\partial s}{\partial x} \\
 +
& x:\quad 4\frac{\partial }{\partial x}\left( \eta \frac{\partial u}{\partial x} \right)+2\frac{\partial }{\partial x}\left( \eta \frac{\partial v}{\partial y} \right)+\frac{\partial }{\partial y}\left[ \eta \left( \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right) \right]+\frac{\partial }{\partial z}\left( \eta \frac{\partial u}{\partial z} \right)=\rho g\frac{\partial s}{\partial x} \\
 +
\end{align}</math>
 +
 +
 +
... with an analogous expression for the ''y'' balance equation. For actually solving the system of equations, it is advantageous to combine all ''v'' gradient terms for the ''x'' balance equation and move them to the right-hand side (and vice versa for the ''y'' balance equation). This gives the final form of the equations to be discretized and solved:
 +
 +
 +
<math>\begin{align}
 +
  & x:\quad 4\frac{\partial }{\partial x}\left( \eta \frac{\partial u}{\partial x} \right)+\frac{\partial }{\partial y}\left( \eta \frac{\partial u}{\partial y} \right)+\frac{\partial }{\partial z}\left( \eta \frac{\partial u}{\partial z} \right)=-2\frac{\partial }{\partial x}\left( \eta \frac{\partial v}{\partial y} \right)-\frac{\partial }{\partial y}\left( \eta \frac{\partial v}{\partial x} \right)+\rho g\frac{\partial s}{\partial x} \\
 +
& y:\quad 4\frac{\partial }{\partial y}\left( \eta \frac{\partial v}{\partial y} \right)+\frac{\partial }{\partial x}\left( \eta \frac{\partial v}{\partial x} \right)+\frac{\partial }{\partial z}\left( \eta \frac{\partial v}{\partial z} \right)=-2\frac{\partial }{\partial y}\left( \eta \frac{\partial u}{\partial x} \right)-\frac{\partial }{\partial x}\left( \eta \frac{\partial u}{\partial y} \right)+\rho g\frac{\partial s}{\partial y} \\
 +
\end{align}</math>
 +
 +
 +
Go to [[Blatter-Pattyn Boundary Conditions]].
 +
 +
Go to [[CISM higher-order dynamics: model solution|solution of the Blatter-Pattyn model]].
 +
 +
Go to [[CISM governing equations and numerical solution]].
 +
 +
[[International Workshop on Ice Flow Modeling, Beijing Normal University (March 21-25, 2011)|Return to main coarse page]]
 +
 +
 +
== References ==
 +
 +
* [http://www.ingentaconnect.com/content/igsoc/jog/2010/00000056/00000197/art00010 Dukowicz, J. K., S.F. Price, and W. H. Lipscomb. 2010. Consistent approximations for ice-sheet dynamics from a principle of least action. ''J. Glaciol.'', '''56'''(197).]
 +
 +
* [http://homepages.vub.ac.be/~fpattyn/papers/Pattyn2003_JGRB.pdf Pattyn, F. A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes, ''J. Geophys. Res.'', '''108'''(B8), 2003].
 +
 +
* [http://qjmam.oxfordjournals.org/content/63/1/73.abstract Schoof, C. and R. Hindmarsh. Thin-film flows with wall slip: an asymptotic analysis of higher order glacier flow models. ''QJMAM'' '''63'''(1), doi:10.1093/qjmam/hbp025].

Latest revision as of 16:09, 14 March 2011

The higher-order dynamics scheme implemented in Glimmer/CISM is based on the so-called "Blatter-Pattyn" model. The starting point for this model is the full Stokes equations


\begin{align}
  & x:\quad \frac{\partial \tau _{xx}}{\partial x}-\frac{\partial P}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=0 \\ 
 & y:\quad \frac{\partial \tau _{yy}}{\partial y}-\frac{\partial P}{\partial y}+\frac{\partial \tau _{xy}}{\partial x}+\frac{\partial \tau _{yz}}{\partial z}=0 \\ 
 & z:\quad \frac{\partial \tau _{zz}}{\partial z}-\frac{\partial P}{\partial z}+\frac{\partial \tau _{zy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial x}=\rho g \\ 
\end{align},


where P is the pressure and τ is the deviatoric stress tensor. The latter is given by


\tau _{ij}=\sigma _{ij}+P\delta _{ij},


where σ is the full stress tensor.


There are a number of ways to argue that, due to the "shallowness" of ice sheets (because the ratio of H/L, where H is the thickness and L is a relevant horizontal length scale, is small), the equations above can be reduced to the following "first-order approximation":


\begin{align}
  & x:\quad \frac{\partial \tau _{xx}}{\partial x}-\frac{\partial P}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=0 \\ 
 & y:\quad \frac{\partial \tau _{yy}}{\partial y}-\frac{\partial P}{\partial y}+\frac{\partial \tau _{xy}}{\partial x}+\frac{\partial \tau _{yz}}{\partial z}=0 \\ 
 & z:\quad \frac{\partial \tau _{zz}}{\partial z}-\frac{\partial P}{\partial z}=\rho g \\ 
\end{align}


The arguments that support this reduction are fairly complex, either based on a variational analysis or an asymptotic analysis. Two of the papers given in the references below (Dukowicz et al.,2010 and Schoof and Hindmarsh,2010, respectively) provide more details on the mathematical background that allows us to state that these equations are "first order accurate" approximations to the Stokes equations.


We continue by noting that the 3rd (vertical) balance equation above can be integrated through the depth of the ice sheet to give an expression for the pressure, P,


P=\rho g\left( s-z \right)+\tau _{zz}(z)


This is simply a statement the the full vertical normal stress is balanced by the hydrostatic pressure (the so-called hydrostatic assumption). This expression can be substituted into the horizontal pressure gradient terms above to remove pressure from the equations. For example, for the x component of velocity we have


\begin{align}
  & x:\quad \frac{\partial \tau _{xx}}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=\frac{\partial P}{\partial x} \\ 
 & x:\quad \frac{\partial \tau _{xx}}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=\frac{\partial }{\partial x}\left[ \rho g\left( s-z \right)+\tau _{zz}(z) \right] \\ 
 & x:\quad \frac{\partial \tau _{xx}}{\partial x}-\frac{\partial \tau _{zz}}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=\rho g\frac{\partial s}{\partial x} \\ 
\end{align}


Now we can use the incompressibility and stress-strain colinearity constraints to rewrite the vertical normal deviatoric stress in terms of horizontal normal deviatoric stresses. This allows us to write the horizontal balance equations in terms of horizontal stress gradients only. Again, for the x direction we have


\begin{align}
  & \tau _{zz}=-\tau _{xx}-\tau _{yy} \\ 
 & -\frac{\partial \tau _{zz}}{\partial x}=-\frac{\partial }{\partial x}\left( -\tau _{xx}-\tau _{yy} \right)=\frac{\partial \tau _{xx}}{\partial x}+\frac{\partial \tau _{yy}}{\partial x} \\ 
 & \frac{\partial \tau _{xx}}{\partial x}-\frac{\partial \tau _{zz}}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=\rho g\frac{\partial s}{\partial x}\quad ...\text{becomes}... \\ 
 & 2\frac{\partial \tau _{xx}}{\partial x}+\frac{\partial \tau _{yy}}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=\rho g\frac{\partial s}{\partial x}\quad  \\ 
\end{align}


Note that at this point we have removed the vertical balance equation entirely; it has been incorporated into the horizontal balance equations using incompressibility. The horizontal balance equations are now:


\begin{align}
  & x:\quad 2\frac{\partial \tau _{xx}}{\partial x}+\frac{\partial \tau _{yy}}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=\rho g\frac{\partial s}{\partial x}\quad  \\ 
 & y:\quad 2\frac{\partial \tau _{yy}}{\partial y}+\frac{\partial \tau _{xx}}{\partial y}+\frac{\partial \tau _{xy}}{\partial x}+\frac{\partial \tau _{yz}}{\partial z}=\rho g\frac{\partial s}{\partial y} \\ 
\end{align}


Next, we want to write these equations in terms of velocities, since that is what we are ultimately solving for. The link between stress gradients and velocities is through the constitutive law for ice (here we'll assume Glen's law) - relating strain rates to stresses - and the definition of the strain rate tensor - relating strain rates to velocity gradients.


\begin{align}
  & 1.\quad \tau _{ij}=B\dot{\varepsilon }_{e}^{\frac{1-n}{n}}\dot{\varepsilon }_{ij},\quad B=B(T) \\ 
 & 2.\quad \dot{\varepsilon }_{ij}=\frac{1}{2}\left( \frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}} \right) \\ 
 & 3.\quad 2\dot{\varepsilon }_{e}=\dot{\varepsilon }_{ij}\dot{\varepsilon }_{ij} \\ 
 & 4.\quad \eta \equiv \frac{1}{2}B\dot{\varepsilon }_{e}^{\frac{1-n}{n}} \\ 
 & 5.\quad \tau _{ij}=2\eta \dot{\varepsilon }_{ij} \\ 
\end{align}


In order, the four expressions above give

  1. Glen's flow law (actually, the inverse flow-law here) where B(T) is a temperature dependent pre-factor
  2. The definition of the strain rate tensor in terms of velocity gradients
  3. The definition of the effective strain rate, \dot{\varepsilon }_{e}, a norm of the strain rate tensor
  4. A definition for the "effective viscosity" (after rearranging some terms in (1))
  5. ... which allows us to write the relationship between stress and strain in a standard "Newtonian" way (but with a non-Newtonian viscosity)


Now taking these definitions into our stress balance equations above and expanding in terms of strain rates, velocity gradients, and effective viscosities, we have (for the x direction only)


\begin{align}
  & x:\quad 2\frac{\partial \tau _{xx}}{\partial x}+\frac{\partial \tau _{yy}}{\partial x}+\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z}=\rho g\frac{\partial s}{\partial x} \\ 
 & x:\quad 2\frac{\partial }{\partial x}\left( 2\eta \dot{\varepsilon }_{xx} \right)+\frac{\partial }{\partial x}\left( 2\eta \dot{\varepsilon }_{yy} \right)+\frac{\partial }{\partial y}\left( 2\eta \dot{\varepsilon }_{xy} \right)+\frac{\partial }{\partial z}\left( 2\eta \dot{\varepsilon }_{xz} \right)=\rho g\frac{\partial s}{\partial x} \\ 
 & x:\quad 4\frac{\partial }{\partial x}\left( \eta \frac{\partial u}{\partial x} \right)+2\frac{\partial }{\partial x}\left( \eta \frac{\partial v}{\partial y} \right)+\frac{\partial }{\partial y}\left[ \eta \left( \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right) \right]+\frac{\partial }{\partial z}\left( \eta \frac{\partial u}{\partial z} \right)=\rho g\frac{\partial s}{\partial x} \\ 
\end{align}


... with an analogous expression for the y balance equation. For actually solving the system of equations, it is advantageous to combine all v gradient terms for the x balance equation and move them to the right-hand side (and vice versa for the y balance equation). This gives the final form of the equations to be discretized and solved:


\begin{align}
  & x:\quad 4\frac{\partial }{\partial x}\left( \eta \frac{\partial u}{\partial x} \right)+\frac{\partial }{\partial y}\left( \eta \frac{\partial u}{\partial y} \right)+\frac{\partial }{\partial z}\left( \eta \frac{\partial u}{\partial z} \right)=-2\frac{\partial }{\partial x}\left( \eta \frac{\partial v}{\partial y} \right)-\frac{\partial }{\partial y}\left( \eta \frac{\partial v}{\partial x} \right)+\rho g\frac{\partial s}{\partial x} \\ 
 & y:\quad 4\frac{\partial }{\partial y}\left( \eta \frac{\partial v}{\partial y} \right)+\frac{\partial }{\partial x}\left( \eta \frac{\partial v}{\partial x} \right)+\frac{\partial }{\partial z}\left( \eta \frac{\partial v}{\partial z} \right)=-2\frac{\partial }{\partial y}\left( \eta \frac{\partial u}{\partial x} \right)-\frac{\partial }{\partial x}\left( \eta \frac{\partial u}{\partial y} \right)+\rho g\frac{\partial s}{\partial y} \\ 
\end{align}


Go to Blatter-Pattyn Boundary Conditions.

Go to solution of the Blatter-Pattyn model.

Go to CISM governing equations and numerical solution.

Return to main coarse page


References