Theory of lattice Boltzmann simulations of glacier flow

Abstract

A lattice Boltzmann technique for modeling Navier–Stokes fluid flow is extended to allow steady-state simulations of glaciers and other slow-flowing solids. The technique is based on a statistical mechanical representation of flowing ice as a set of particles (populations) which translate and collide on a face-centered cubic lattice. The average trajectories of the populations give the velocities of the ice at any point in the glacier. The method has considerable advantages over other techniques, including its ability to handle complex realistic geometries without additional complications to the code Examples are presented for two-dimensional simulations.

Introduction

A common, if unstated, goal of many glaciological field programs is to collect surface-velocity and geometry information suitable for numerical inversions for conditions at the unobserved bed. Two classic examples from the recent past (among many) are Columbia and Variegated Glaciers with volumes of published information and numerous papers devoted to numerical inversions (e.g. Balise and Raymond, 1985; Meier and others, 1985; Raymond and others, 1987; Van der Veen and Whillans, 1993). In fact, during the last century, dozens of other glaciers and ice streams have also been repeatedly to gather information suitable for inversions aimed at elucidating sliding laws and mechanical controls at or near the bed. Typical numerical procedures have focused on finite-element and finite-difference formulations of the inversion problem, and with the advent of work stations, these methods have become increasingly accurate, especially in their descriptions of the deformational components of flow. However, the model’s efficiency and accuracy in describing complex geometries have lagged far behind the more impressive developments in continuum flow. Finite-element techniques can take advantage of new methods for automatic generation of appropriate modeling meshes and finite-difference techniques can utilize specialized curvilinear coordinate systems. In both cases, however, the increased ability to deal with realistic geometries has forced a corresponding increase in computational and coding complexity.

Within the last few years a number of exciting advances in statistical mechanics have opened up new approaches to modeling low-viscosity fluid flows around complex geometric structures. One of the most recent and successful of these new techniques is the lattice gas automaton (LGA) which excels in its ability to handle complex boundary conditions with a surprisingly simple code (Frisch and others, 1986). LGA models represent a fluid as a set of colliding gas particles which on average reproduce incompressible Navier-Stokes fluid flow. Boundary conditions are handled by simple particle-collision rules which reproduce the desired fluid velocities at any point. Using LGA, for example, flows through complex, even fractal, porous media geometries have been successfully modeled at a level of detail not previously possible (e.g. Chen and others, 1991). Likewise, simulauons of flow around dissolving miscible sohds and simulations of turbulent flow around arbitrary obstacles has become routine (e.g. d’Humières and Lallamand, 1987; Lawniczak and others, 1991).

These automaton techniques have been slow to excite interest in disciplines which focus on the slow-creeping flow of solids, like glaciers, primarily because the bulk of LGA research has developed models suitable for low viscosities and high Reynolds number simulations. These constraints, however, can be removed, and using LGA to simulate glaciers requires only a shift in modeling philosophy. Rather than focusing on finite-difference or finite-element models of the macroscopic equations of state, LGAs focus on building much simpler models formulated at a microscopic level. The microscopic models are constructed to conserve mass, momentum and energy but do not attempt completely faithful reproductions of the microscopic world. This is a significant departure from traditional modeling philosophy but, although inexact, this type of simulation is successful because hydrodynamics and thermodynamics describe large-scale systems which behave independently of thar precise nucroscopic formulation (Salem and Wolfram, 1986). So, while the macroscopic models are typically limited to a small number of spatial “cells” where information on the flow state is updated based on the complex continuum equations, these microscopic models instead use tens of thousands of cells which are updated by simple rules chosen for their computational efficiency. Proper continuum behaviour is reconstructed from large-scale averages of the many microscopic cells. In other words, very complex systems are studied by looking at simplified microscopic models which reproduce the correct mathematical behaviour at macroscopic scales.

The shift from a macroscopic to a simplified microscopic paradigm will allow significant advances in our ability to model the flow of glaciers, ice sheets and other solids around complex realistic obstacles. In the following sections, therefore, we describe a variant of LGA, called the lattice Boltzmann method, and demonstrate a new technique for implementing constitutive relationships appropriate for ice and other non-Newtonian slow-flowing solids. The lattice Boltzmann method described here is appropriate for isothermal steady-state problems or for time-varying velocities with fixed glacier geometries. The technique can he used to invert surface velocities for information at the bed or to extract velocity and stress profiles from complicated, realistic glacier channels. This approach represents the first step towards a more general lattice Boltzmann model of time-dependent and thermally inhomogeneous glacier flow.

The Lattice Boltzmann Method

Historically, the lattice Boltzmann method evolved from LGA and is most easily understood in the context of its automaton predecessors. Automatons correspond closely with von Neumann’s original concept of a computer and are any machine with a finite set of input and output states and a fixed finite set of rules which map each input to an output (Lewis and Papadimiviou, 1981, p. 222). A simple example uses finite-length binary strings as input and output states; a look-up table converts a binary input string to a binary output string. The set of rules, or the look-up table, is the automaton’s program.

Lattice gas automatons have input and output states which correspond to particle positions and velocities on a triangular lattice. At each site on the lattice, up to six particles of unit mass are assigned unit velocities e _i in one of the six lattice directions (i = 1, …, 6) (Fig. 1). At each point on the lattice there is at most one particle with a given velocity (in other words, no two particles are allowed to have the same momentum at the same lattice site). If n_i = 1 and n_i = 0 represent occupied and unoccupied momentum states, then the particle positions and velocities at lattice site x and time t are given by the binary string n(x, t) = (n ₁ n ₂ n ₃ n ₄ n ₅ n ₆). Figure 1 illustrates n(x, t)= (100001).

Fig. 1. Momentum state n(x,t)= (100001) on a triangular lattice.

At each time step in the model’s evolution, particles in the LGA are translated one lattice unit in the direction of their velocity. Particles which arrive at the same site are said to collide and are redistributed according to a fixed set of rules which conserve mass and momentum (these are the automaton rules which convert the binary input to the binary output states). After colliding, the particles are translated again, and the process repeats. Figure 2 illustrates the translation and collision rules for two and three particles.

Fig. 2. An example of the translation and collision process for LGA particles. The three-way collision at site x re-arranges velocities in the only manner which will conserve momentum and mass. The two particles collidingcollinearly at site y have zero total mometum and can be redistributed in any other collineor directions (the choice can be made randomly). All other possible collision configurations are ignored and the velocities are left unchanged (for example, grid site z)

This basic lattice gas automaton was first proposed by Frisch and others (1986) in a paper in which they demonstrated that the Navier-Stokes equation is reproduced in the continuum limit by the movement of the particles on the triangular lattice. If the local density ρ and local velocity v of the fluid at any point in time are given by

(1)

and

(2)

then they showed that macroscopic quantities like the fluid velocity are obtained by spatial averages of v. Boundary conditions are implemented by adjusting the directions of any particles which collide with a boundary; and at each time step, body forces, like gravity, can be specified by altering the direction (momentum) of random particles to align with the force (Kadanoff and others, 1987). A number of review papers have given details of the derivations which lead from the collision rules to the full Navier Stokes equations of fluid flow (e.g. Wolfram, 1986; Frisch and others, 1987).

LGA models of fluid flow are stable because the lookup tables are exact Boolean operators (round-off errors arc eliminated). The translation and collision process is also local in nature, so that the updates can he performed simultaneously at all the nodes, making LGA well suited for fast parallel processing. The discrete nature of the particle momenta, however, adds noise to the system, requiring large spatial averages of the particle trajectories to get reasonable fluid approximations. This noise clearly outweighed the advantages of an LGA implementation of the glacier-flow problem constructed by the authors, and along with a problem in adjusting the viscosity (controlled by the selection of specific collision rules) motivated others to develop the original lattice Boltzmann technique (Higuera and Jiménez, 1989).

The lattice Boltzmann method considers populations of particles (0 ≤ n(x, t) ≤ 1) rather than individual particles on the lattice. In essence, this is the same as simultaneously running a large ensemble of LGA simulations and then determining the probability of a particle having a particular momentum state (i.e. probability of occupying a particular lattice site with a particular velocity). (In statistical mechanics, this collection of LGA simulations would be referred to as the canonical ensemble which has a Boltzmann distribution of possible states, hence the name.) Practically speaking, the lattice Boltzmann technique averages out the noise in the system while maintaining the fast parallel architecture. However, the arithmetic is now real, rather than Boolean, and the ultimate stability of the LGA is compromised, though stable behavior is still possible for almost all flow scenarios of interest. (Benzi and others (1992) discussed details of stability to numerical dispersion, discretization and other types of error.)

The lattice Boltemann collision rules can be implemented in the same manner as the LGA rules, although a more efficient method is discussed here. In the original LGA, there were six unit velocities e _i (i = 1, …, 6) on a triangular lattice, but more recent analyses (Skordos, 1993) have shown that the Navier-Stokes equations can he reproduced with two velocities on a more easily coded square lattice (velocities along edges and diagonals). In particular, for i = 1, …, 8

(3)

In addition, e ₀ = (0, 0) corresponds to a “rest” particle with zero velocity. Also, rather than using look-up tables, the collision step can he more compactly represented as a relaxation from an equilibrium state. Let 0 ≤ N_i (x, t) ≤ 1 be the probability of finding a particle at node x and time t with unit velocity e _i and let N_i ^eq be the equilibrium distributions. Then, for a unit time step (Δt = 1),

(4)

The lefthand side of Equation (4) is just the rate of change of the particle distributions due to collisions. τ is a unitless relaxation parameter which controls this change and adjusts the rate at which the equilibrium is approached.

The equilibrium states N_i ^eq are derived from a Chapman-Enskog expansion of Equations (1), (2) and (4) (Skordos, 1993) and for a square lattice are given by

(5)

The density ρ and macroscopic velocity v are given by Equations (1) and (2) without averaging.

The relaxation time parameter is also derived from the Chapman-Enskog expansion and is related to the kinematic viscosity v (measured in lattice units squared per time step) by

(6)

τ has a lower bound of ½ to prevent unphysical negative viscosities. Intuitively, the relationship between τ and v is a consequence of the relaxation parameter’s influence on the rate of change of particle distributions at each collision step Thinking of the Boltzmann method as an ensemble of LGAs, a larger τ decreases the change in Ν with time (Equation (4))or, in other words, decreases the probability of a collision in the ensemble of LGAs. When collisions are less probable, the mean free path of a particle is increased and Maxwell (1859) demonstrated that the mean free path of particles in a gas is proportional to the viscosity (Huang, 1987). Therefore, increasing τ increases the mean free path and increases the viscosity.

Non-Newtonian Flow

The ability to adjust the viscosity with the relaxation parameter is our basis for constructing non-Newtonian suitable for modeling glacier flow. When τ has a constant value, the viscosity is fixed and stress is linearly related to strain. In other words where and are components of the deviatoric stress and strain-rate tensors and μ = ρν is the dynamic viscosity (e.g. Mase, 1970). Water is the classic example of such a linear or Newtonian fluid Ice, magma, metals and many other solids are characterized, however, by an effective viscosity, μ _eff, which is a function of the strain rate (in addition to other dependencies such as chemical impurities and temperature). In fact, for ice, magma and most metals, the effective viscosity is a power of the second strain-rate invariant,

(Turcotteand Schuhert, 1982; Paterson, 1994), i.e.

(7)

where

(8)

for some positive n.

To handle these generalized viscosities, τ is simply expressed as a function of the effective viscosity,

(9)

Now the strain rates at each lattice site can be estimated differences of the velocities at their neighbors. This determines the effective dynamic viscosity μ _eff and an appropriate τ at each node on the lattice. Aharonov and Rothman (1993) have developed a similar but more complicated scheme for non-Newtonian flow based on an earlier implementation of the lattice Boltzmann method. Our results compare favorably with their work and require less complicated coding.

Eliminating Compressibility and Inertial Effects

Solids typically have kinematic viscosities on the order of 1 × 10¹⁵ m² s⁻¹ about 21 orders of magnitude higher than the typical viscosities for a fluid (1 × 10 ⁻⁶ m⁻²s⁻¹). For a lattice Boltzmann model, this has two important consequences for flow simulations of a solid. First, the mean free path of a particle is so long that infeasibly large model grids are necessary for realistic simulations. This can only be countered by taking advantage of dynamic similarity and using dynamically equivalent flow scenarios with much smaller viscosities. Secondly, large viscosities imply negligibly small inertial terms (i e for Newtonian flow, the v. Δv term in the Navier-Stokes equation should be set to zero). However, these two consequences tend to work against each other; decreasing the viscosity causes increasingly dominant inertial effects in the model, leading to high Reynolds numbers and turbulent flow. To prevent this behavior, it is best to eliminate the intertial terms from the glacier-model formulation.

To do this, consider the macroscopic flow behavior of the lattice Boltzmann model as derived using the Chapman-Enskog procedure, which is essentially a Taylor expansion of Equation (4). The expanded N_i are rewritten in terms of ρ and ρ v using the mass and momentum relationships given in Equations (1) and (2). The equilibrium distribution is treated as an expansion about small velocities (Hou, 1995),

(10)

with higher-order terms neglected. The general parameters A_i, B_i, C_i. and D_i are then selected to ensure that to first order the expansion gives the Navier-Stokes equations (see Equation (5)).In our case, we choose a different set of parameters which gives the same solution but eliminates the inertial terms. A_i and B_i remain the same as before, but C_i = D_i = 0. Intuitively, the terms in Equation (10) which are non-linear in v have been set to zero because the Navier-Stokespdes are linear in v except for the inertial term, v. Δv.

In our model, one additional correction is made to ensure incompressibility Without some further manipulation, the equilibrium distributions given in Equations (5) and (10) predict that Δ.(ρ v) = 0 for steady-state flows with velocity v. Truly incompressible flows require Δ . v = 0. Therefore, Hou (1995) has suggested redefining the velocity to be

(11)

Using the same Chapman-Enskog expansion as outlined above, this leads to a correction in the equilibrium given by

(12)

A more complete motivation of the incompressibility corrections has been explained in Hou (1995).

Generalization to Three Dimensions

A three-dimensional generalization of the lattice Boltzmann model is straightforward. Instead of a square lattice, the particle populations move on a face-centered cubic lattice with six nearest neighbors at a distance of |e _i ∆t| = 1 lattice unit and eight next-nearest at a distance of All other equations are unchanged except for the equilibrium distributions which become

(13)

for i = 0 to 14.

Scaling Analysis and Non-Dimensionalization

In order to run the lattice Boltzmann model at low viscosities and have dynamically similar behavior to a highly viscous glacier, a set of non-dimensional numbers needs to be derived from the macroscopic equations of glacier flow. Τo guarantee similar dynamics, these nondimensional numbers must be the same for both the high- and low-viscosity scenarios. For Newtonian flows, the appropriate non-dimensional numbers are derived by dividing each term of the Navier-Stokes equation by one of the other terms (Welty and others, 1984). Reynolds number, for instance, is the ratio of the inertial term to the viscous term However, the equations of glacier flow are complicated by a non-linear rheology (Equation (7)) and cannot be written as a single partial differential equation. In that case, ratios of terms are impossible and a more sophisticated scaling analysis is necessary to derive the non-dimensional parameters (Logan, 1987).

Scaling analyses start with a set of relevant variables {V ₁, V ₂, ... V_i } and replace each by a rescaled quantity The rescaled variables are substituted into the relevant equations (describing glacier flow in our case) and the scaling parameter λ is factored out. This factoring ensures that the equations are unchanged by the rescaling and is used to derive a set of constraints on {m ₁, m ₂, ... m_i }. These constraints determine how each variable must be related to its rescaled quantity and therefore determine how different scale simulations must he related to each other in order for the relevant equations to describe the same physical process. In other words, the constraints provide the relationships necessary to construct non-dimensional quantities.

Consider the equations describing isothermal steady-state glacier flow in two dimensions. From the definition of strain rates,

where x and z are along coordinate directions perpendicular and parallel to gravity, and u and v are deformational velocity components in the x and z directions. From these definitions,

(14)

which is a compatibility requirement ensuring that the velocities are unique (e.g. Mase, 1970, p. 92). The strain rates are related to deviatoric stresses by Glen’s flow law for ice (Equation (7)), which can be written explicitly as

(15)

and

(16)

with and . (σ_ij is a stress-tensor component.) Conservation of momentum requires

(17)

and

(18)

where g is gravity and α is the angle of the ice surface relative to the horizontal.

The variables in these equations can be rescaled by the following substitutions, x = λ ^a x, z = λ ^b z, σ _xx = λ ^c σ_xx , σ _xz = λ ^d σ_xz , σ _zz = λ ^e σ_zz , A = λ ⁱ A, ρ = λ ^j ρ, g = λ ^k g, u = λ ^p u and υ = λ ^q v, where the overbar represents a scaled quantity and λ is an arbitrary scaling factor. After substitution, the rescaled equations are consistent with the original equations only if the As can be factored out. This places a set of constraints on the scaling exponents. For example, using the chain rule, (∂σ _{x x} )/(∂ x ) = λ ^c–a (∂σ _{x x} )/(∂x) and (∂ σ _{x z} )/(∂ z ) = λ ^d–b (∂σ _{x z} )/(∂z), so

(19)

This gives back the original Equation (17) if and only if λ ^c–a = λ ^d–b = λ ^j+k or c – a = d – b = j +k.

Similar constraints derived from each of the other flow equations leads to a set of relationships which can be solved to give a = b, c = d = e = a + j + k and p = q = a + i +n(a + j+ k). Using λ ^a = x /x etc., these solutions imply that

(20)

(21)

and

(22)

Equations (20) through (22) give constraints that will be satisfied when two simulations are dynamically similar.

By solving for the original unsealed variables, Equations (20) through (22) can be rewritten as non-dimensional number

(23)

(24)

and

(25)

Each of the variables now represents characteristic values (e.g. the glacier thickness for x). Note that the characteristic shear stress is given by the basal shear, or ρgx sin α, so the non-dimensional parameter π ₂ just requires the slopes to be the same in dynamically similar simulations π ₃ requires geometric similarity and π ₁ is a generalization of Reynolds number divided by Froude’s from Newtonian flow (n = 1). Note that the inclusion of a third dimension in the glacier-flow equations does not alter π₁ or π ₂ and leads to the obvious generalizations cf geometric similarity. The axes orientations are also irrelevant to the derivation

As a typical example, fix n = 3 and let the glacier thickness be x = 100 m, let the characteristic velocity be u = 1/100000 m s⁻¹ g = 9.8 m s⁻², ρ = 900 kg m⁻³ and A = 1 × 10⁻²⁴ m³ s⁵ kg⁻³. Then π₁ , = 6.86. A dynamically similar model could have the values x = 30, u =0.1, g = 4.73 × 10⁻⁴, ρ = 2 and A = 1000 where each variable is measured in model units. These values are typical for the lattice Boltzmann model and note that the viscosity is kept relatively small, as required.

Two-Dimensional Glacier-Flow Simulations

A glacier-flow simulation starts by occupying each lattice site and lattice direction with some initial population of These populations are then translated in the direction of their velocities. Populations arriving at the same lattice site are redistributed by the rules which conserve mass and momentum (Equations (4) and (12)). As in the lattice gas models, gravity is implemented by changing the momentum of particle populations by some small increment at each lattice site at each time step The translations, redistributions and body-force momentum ddjustments iterate until the system reaches a steady state. In our simulations, a change in velocity (in a least-squares sense) of less than 1 × 10⁻⁵ between iterations is required for steady state. Velocities are calculated from Equation (11).

Equation (12) gives the basis for specifying velocity-boundary conditions (Skordos, 1993). Rather than calculating the velocity from Equation (11), the velocity at any boundary site is fixed to the desired value (for example, to that of a surveyed velocity). At each iteration, the specified velocity is used to calculate an equilibrium distribution via Equation (12). The specified equilibrium distribution is then used to adjust the populations in the collision step (Equation (4)).Specifying zero velocity, for example, leads to a no-slip condition. A free-slip (or frictionless) condition is also possible by adjusting Equation (4) so that the angle of incidence equals the angle of reflectance for any populations approaching the free-slip boundary. Note that both specified velocities and free-slip conditions can be used at different sites in the same simulation.

Two examples of this simulation process are presented below. An implementation of the three-dimensional model is in progress and examples will be presented in future publications.

Laminar flow

The simplest possible glacier is a uniform thickness of ice sitting on an inclined plane. Assuming a no-slip condition at the bed and a free-slip condition at the surface, the velocity profile is given by the laminar-flow solution, υ = 0 where u_s is the surface velocity and is z measured perpendicular to the surface, positive downwards (Paterson, 1994, p. 251). For n = 1, A = 3, ρ = 1.0, g = 2.8 × 10⁻⁵. α = 0.46 rad and a bed at z = 30 this gives u(z) = 0 034 - 3.7 × 10⁻⁵ z ². For this case, π₁ = 2.22, so one choice of corresponding physical values is A = 1 × 10⁻³ m s² kg⁻¹, ρ = 900 kg m⁻³, g = 9.8 ms², a depth of 1000 m and a slow surface velocity of 0.0125 m a⁻¹ Figure 3 compares the results of the Boltzmann model to the analytical expression. A similar simulation with n = 3, g= 1.2 × 10⁻³ and A = 1000 illustrates the characteristic flattening of the velocity profile caused by a strain-thinning solid (Fig. 3). In this case, π₁ = 22.86, so one choice of corresponding physical values is A = 1 × 10⁻²⁴ m³ s⁵ kg⁻¹, ρ = 900 kg m⁻³, g = 9.8 m s⁻². a depth of 100 m and a surface velocity of 94.66 m a⁻¹.

Fig. 3. Laminar flow with n — l (top) and n = 3 (bottom). Point data are from the lattice Boltzmann simulations and solid lines are the analytical solutions.

Valley-glacier flow

Figure 4 shows non-Newtonian (n = 3) flow through a more typical two-dimensional cross-section of a valley glacier. A no-slip condition is specified at the bed and surface velocities are fixed at the surface (roughly downwards in an accumulation zone and roughly upwards in an ablation zone). Other parameters (and putential choices for characteristic physical values) are the same as above. Notice that the bed is undulating and overdeepened near the terminus and that the cross-section has an icefall and a steep bedrock cliff near the equilibrium line. These geometries can cause significant meshing difficulties for non-Newtonian finite-difference and finite-element techniques but pose no difficulty for the lattice Boltzmann model.

Fig. 4. Velocity vectors calculatedby the lattice Boltzmann model for a two-dimensional vertical cross-seclion through an irregularly shaped valley glacier. Gravity is acting vertically downwards. The small dots are lattice sites which are not a part of the glacier and have zero velocity. The vectors are magnified 100 times and are measured in lattice units per time step. The grid is 50 by 30 lattice units.

Conclusions

The code for a lattice Baltzmaim model is simple and short, with an intuitive motivation — the translation and collision of particles. Glacier-flow geometries can be as complicated as necessary with no additional complexity imposed on the model code. The same technique can be used to model the slow creeping flow of any solid, such as glass, cold tar, salt diapers or lava. The only restrictions are a quantifiable constitutive relationship (substituted in place of Equation (8)), flows which have negligible inertia, steady-state behavior and isothermal conditions.

The preliminary success of the lattice Boltzmann method for steady-stale isothermal glacier flows suggests that it would be profitable to focus additional research on fully time-dependent three-dimensional thermomechanical Boltzmann techniques. The current model can already handle time-varying boundary conditions, so the primary difficulty is allowing free-surface geometries which can evolve with time. The authors are currently developing this free-surface condition as well as techniques for incorporating the thermodynamic equations of state.