One of the recent blog posts discussed the OGGM time stepping scheme. In short, it explains that OGGM still uses an ad-hoc CFL condition to determine the time step. An “optimal” value, currently in use, was determined through many empirical tests, making a compromise between speed and error growth. However, since this value was not determined by a rigorous mathematical treatment:
- it does not ensure numerical stability for all glaciers: some would need a lower CFL condition to prevent error growth,
- it is not optimal: for some glaciers, we could possibly use a higher CFL value and thus reduce the computation time.
So why don’t we just crunch the numbers and determine the optimal (rigorously speaking) CFL value? Well, it’s because stability analysis of our time stepping scheme is not that easy. In this blog post we will look at why that is.
Some math & physics basics
As so many other simple glaciological models, OGGM is also based on the shallow ice approximation. In this formulation, the building blocks of glaciers are vertical columns of ice, each moving with a velocity dependent on the ice thickness and the inclination (gradient) of the ice surface (a bit more detail can be found e.g. in a textbook by Hutter, K. (1983).
To model the movement of the ice, we have some flexibility in how we formulate the governing equations. Shallow ice models usually solve diffusion-type equations - this is because generally, these cause fewer numerical issues than advection-type equations.
A short reminder: in the diffusion equation, the time derivative of a variable $\psi$ depends on the divergence of its gradient. On the other hand, in the advection equation, the time derivative of a variable depends on its flux divergence:
where $\psi$ is the variable of interest, $D$ is a diffusion coefficient and $u$ is velocity of the flow. In diffusion-based models $\psi$ is the surface elevation.
However, a diffusion formulation does not provide enough flexibility for OGGM; we want to allow for changing widths of the glaciers and various bed shapes (rectangular, trapezoidal or parabolic). Thus, in the OGGM model we choose to solve an advection-type equation, which gives us a bit more flexibility in this regard. In OGGM, the advected variable is the glacier cross-section (denoted by uppercase $S$):
$\frac{\partial S}{\partial t} = w\dot{m} - \nabla \cdot \left(uS \right)$
where $w$ is the width of the glacier at its surface, $\dot{m}$ is the mass balance and $u$ is the column integrated velocity. In short, the model advects the glacier cross section (and thus mass) along the flow. This formulation allows us to solve the same equation for each bed shape. More details about this formulation can be found in the OGGM documentation.
Since other models use different formulations and solve other equations, we can’t just save ourselves the work and look up the stability conditions. And so we will have to try and tackle this by ourselves.
The usual approach: von Neumann stability analysis
The most common and basic method for determining the CFL criterion is the so-called von Neumann analysis. This approach requires first decomposing the variable of interest into a discrete Fourier series, and then mathematically determining the conditions under which the amplitude of each separate Fourier mode does not grow over time. However, this method can be used only under certain conditions ([Wesseling, P. (2009)], Sec. 5.8).
The first limitation is that this method works only for linear equations with constant coefficients (in our case, that would mean that $u$ is constant in space); in non-linear equations, there is an interaction between various Fourier modes, and thus the modes cannot be each treated separately. Second, again because of the Fourier decomposition, we require either a periodic boundary condition, or an infinite spatial domain. The third requirement is constant grid spacing. Let’s look at all these conditions one by one.
Linearity
In the shallow ice equation, the velocity $u$ depends among others on the ice thickness $h$ and the surface inclination $\alpha$ (defined as $\alpha = -\frac{\partial s}{\partial x} = -\frac{\partial}{\partial x}(h+b)$, where lowercase $s$ is the ice surface elevation, and $b$ is the bedrock elevation). Thus, the velocity is variable in space:
$ u = f_d (\alpha \rho g)^{n} h^{n+1} + f_s (\alpha \rho gh)^{n} h^{n-1} . $
(Again, for more details and definitions, go here.)
Furthermore, for all bed shapes, the cross section $S$ is a function of the ice thickness $h$:
$\qquad S = w_0h + \frac{\lambda h^2}{2}$ (trapezoidal)
$\qquad S = \frac{2}{3}wh=\frac{4h^{3/2}}{3P^{1/2}} $ (parabolic)
Here, the width $w$ of a rectangular bed is a given constant. For the trapezoidal bed, $w_0$ is the bed width at the valley floor and $\lambda$ determines the wall angle. The bed shape of a parabolic glacier is further determined by a shape parameter $P_s = \frac{4h}{w^2}$, so that the cross section is again determined only by one constant parameter and the ice thickness.
For each bed shape we can invert this relationship between the cross-section and the ice thickness, and say that $h=h(S)$:
$\qquad h = \frac{\sqrt{w_0^2 + 2 \lambda S}-w_0}{\lambda}$ (trapezoidal)
$\qquad h = \left( \frac{3}{4}P_s^{1/2}S \right)^{2/3}$ (parabolic)
Because of the linearity condition, it is important to stress again that cross-section and ice thickness are not independent for any bed shape. Then, when the shallow ice equation is written out in full, it is easy to see that the term $\nabla \cdot (uS) $ is indeed pretty complicated and non-linear because of the dependence of $u$ on $\alpha$ and $h$ and thus also indirectly on $S$:
$ \frac{\partial S}{\partial t} = w\dot{m} - \nabla \cdot \left[ \left(f_d (\alpha \rho g)^{n} h^{n+1} + f_s (\alpha \rho gh)^{n} h^{n-1} \right) S \right] $
For now, basal slipping is ignored in the model (i.e. $f_s = 0 $), but the first term still remains non-linear: the linearity condition is not fulfilled. This fact is based purely on the differential equation we’re solving and has nothing to do with its numerical implementation.
Boundary conditions
This condition for the von Neumann analysis is also not fulfilled in our model - we do not use periodic boundary conditions. At the uppermost point of the glacier, we do not allow flow into the numerical domain.
Grid spacing
This requirement is the only one fulfilled - the grid spacing indeed does not change over the spatial domain in the OGGM numerical model.
Thus, even not taking the numerical implementation into account, the basic von Neumann method is not usable for determining the CFL condition for our numerical scheme.
Perturbation method
The next thing to try would be to perform the stability analysis using the perturbation approach.
This method is nicely illustrated by [Budd, W. F., and Jenssen, D. (1975)], who apply it to an advection equation similar to ours (albeit less complicated). The central assumption of this method is that the glacier is in steady state - to analyze stability, we then introduce a small perturbation in the variable of interest (ice thickness in the case of Budd and Jenssen, or the glacier cross-section in the case of OGGM) and see under what circumstances this perturbation grows in time, and under what circumstances it decays.
The first step is to separate the variable of interest into a steady state part and a perturbation, and to determine the linearized differential equation for the perturbation growth in time. This step is still independent of the numerical scheme. However, we need to consider the interplay between various variables - a perturbation in the cross-section $S$ will of course lead to perturbations in both the ice thickness $h$ and the surface gradient $\alpha$. Thus, we have to express both $h$ and $\alpha$ in terms of the cross-section. Because of the different relationships between $S$ and $h$ for various bed shapes (see the section above), already at this stage of analysis we can see that the linearized equation and thus also the CFL condition will vary between the three cases.
Once we have the linearized differential equations for error growth, the numerical scheme comes into play and von Neumann analysis can be used to determine the stability condition. In principle, this is doable - however, even for a simpler case considered by Budd and Jenssen, the math is already quite complicated, and the resulting expression for a stability condition needs further assumptions and simplifications to be practically usable. With our equation, the analysis would be very challenging for the parabolic and trapezoidal bed shapes.
Furthermore, even if we did find the CFL condition for our equation and the numerical scheme with this approach, we would have to keep in mind the central assumption (i.e. that the glacier is in steady state), which is rarely fulfilled.
Other ways to perform stability analysis
Or, some of the dead ends I found on the way
In literature on numerical modeling you can find a few more ways of performing stability analysis than just the ones mentioned above. However, whether we can actually use the methods is restricted by the equation we’re solving as well as by the numerical implementation of it. For completeness, I would like to briefly mention some of the methods that I stumbled upon, and explain how they work and why we cannot apply them to our problem.
Local von Neumann analysis
A very well written textbook by [Hirsch, C. (2007)] discusses this method under the section “Stability analysis for non-linear problems” (Ch. 8).
When the coefficients are variable ($u \neq$ const.), this method can be used to determine the necessary (but not sufficient) conditions for stability. This method freezes the coefficient in space, and determines the stability conditions locally. The most stringent requirement are then chosen and applied to the entire domain.
However, this method could be used only if we solved an equation of the form
$ \frac{\partial S}{\partial t} = u(x)\frac{\partial S}{\partial x} $
where $u$ is variable is space, but it is independent of $S$. However, in our case we have an equation of the form
$ \frac{\partial S}{\partial t} = \frac{\partial \left( u(S)S \right)}{\partial x}, $
and so we cannot apply this method. As Hirsch states, this method is only applicable “for linear problems with non-constant coefficients”.
Matrix/eigenvalue method
I read about this method before I realized that I have to look for methods that work with non-linear equations - I found it relevant because it can be used for problems with non-periodic boundary conditions. And although we can’t use this method, I decided to mention it here because it seems to be pretty widely used. (If we wanted to use this method, we would have to linearize our equation first again.)
As the name already implies, the method relies on linear algebra to find stability conditions. We start with the linear initial boundary value problem $ \frac{\partial u}{\partial t} = L(u), $ where $u$ is some variable of interest, and $L$ is a space differential operator acting on $u$. This equation can be discretized and reformulated as a system of differential equations, one for each grid point: $\frac{dU}{dt} = SU + Q. $ Here, $U$ is the vector of $u$-values at grid points, and $S$ is a space discretization matrix - it determines how these values are combined to form spatial derivatives at each point. Finally, $Q$ contains the non-homogeneous terms and boundary values. The structure of $S$ depends on the chosen numerical scheme.
The stability analysis is then based on the eigenvalue structure of the matrix S - if all the eigenvalues don’t fall below a certain value, the scheme is unstable. However, finding the eigenvalues in the first place can be pretty challenging. A good detailed explanation can again be found in Hirsch (Ch. 10).
Take home points
- Stability analysis for non-linear problems with non-periodic boundary conditions is hard.
- The only usable analysis approach that I found is the perturbation method - however, until we do this, the CFL will stay ad hoc.
References
- Hutter, K. (1983). Theoretical Glaciology, Mathematical Approaches to Geophysics, D. D. Reidel.
- Wesseling, P. (2009). Principles of computational fluid dynamics (Vol. 29). Springer Science & Business Media.
- Budd, W. F., & Jenssen, D. (1975). Numerical modelling of glacier systems. IAHS publ, 104, 257-291.
- Hirsch, C. (2007). Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics. Elsevier.