1.1. Background

Modern studies of the behavior of glaciers, ice sheets and ice streams rely heavily on both observations and physically based models. Data acquired via remote sensing provide critical information on geometry and movement of ice over large sections of the Antarctic and Greenland. Although these datasets represent significant advancements in terms of spatial coverage and the variety of processes we can observe, the physical systems to be modeled are nevertheless imperfectly observed. Uncertainties associated with measurement errors are present and physical models are also subject to uncertainties. There is therefore a need for combining observations and models in a fashion that incorporates uncertainty and quantifies its impact on conclusions.

The goal of combining models and observations is hardly new in glaciology (e.g. Reference MacAyealMacAyeal, 1993; Reference Arthern and HindmarshArthern and Hindmarsh, 2003; Reference Gudmundsson and KnightGudmundsson, 2006) nor in the broad areas of the geosciences (e.g. data assimilation as practiced in numerical weather forecasting). We address the goal by formally modeling the uncertainties present, then using Bayes’ theorem to deduce information about all unknowns in the data. We focus on the development of statistical models with strong reliance on physical modeling, a strategy Reference BerlinerBerliner (2003) called ‘physical–statistical modeling’. This is different from traditional physical modeling, with perhaps data-based parameter estimates, and traditional statistical modeling which possibly relies upon qualitative physical reasoning.

An in-depth illustration based on data from the ‘Northeast Greenland Ice Stream’ (NEGIS) is presented. The physical models used are simplified approximations, but are familiar and provide an accessible arena for exemplifying the physical–statistical approach. Specifically, we develop statistically enhanced versions of a simple model of basal shear stress as the sum of driving stress and a corrector process representing other effects, such as lateral resistive stress. Inference for these unobservable stresses is a primary goal here. Bayesian inference is carried out using recently acquired datasets of ice thickness, surface elevation and surface velocity.

The introduction and uncertainty analysis of the corrector process is a crucial feature of our approach. Rather than introducing the correctors, it is natural to suggest that the effects they represent should be formally modeled. This would require more assumptions and/or more data, thereby introducing additional sources of uncertainty. In any case, we do not expect any implementable physical model to remove concerns regarding its applicability nor its uncertainty. Hence, the ability not only to combine models and observations, as discussed by Reference MacAyealMacAyeal (1989), Reference BlatterBlatter (1995) and Reference Joughin, Fahnestock, MacAyeal, Bamber and GogineniJoughin and others (2001), but also to analyze indicators of model validity and uncertainty is paramount.

The primary output of a Bayesian analysis is a posterior distribution, namely the joint probability distribution for unknown quantities conditional on the observed data. Even in our simple illustration, explicit presentation of their joint distribution is not feasible. Hence, a key aspect of Bayesian analysis in such settings is the ability to generate realizations or ensembles from the posterior distribution; the posterior is then studied through statistical summaries of such ensembles. A key by-product of such techniques is the ability to provide ensembles of initial conditions for use in a dynamic forecasting model. That is, one application of the modeling here is as a method for combining models and observations in support of full dynamical prediction or forecasting.

1.2. Glaciological motivations

That glaciers flow under the force of gravity means that important factors in determining velocities include the constraints of the constitutive relationship (Reference Whillans, Van der Veen and OerlemansWhillans, 1987) and ice thickness, in combination with resistive forces acting along the sides and at the base of the glacier. The lithostatic stress at some depth in the ice equals the weight of the ice above that depth, and horizontal gradients in this stress drive ice-stream flow. The flow of the ice stream is impeded by resistive stresses given by the difference between the full stress and the lithostatic component. A separation is therefore made between action or gravitational forces and reaction or resistive forces. After partitioning the full stress into lithostatic and resistive components, substitution into the horizontal balance equation and integration over the full ice thickness gives the following equation for driving stress along the flow (e.g. Reference Van der VeenVan der Veen, 1999):

(1)

where τ _bx is basal shear stress and τ _Rx represents resistive stresses arising from gradients in longitudinal stress and lateral drag. Defining the driving stress as

(2)

where s is ice-surface elevation, H is the ice thickness, ρ is the density of ice and g is the gravity constant, we have

(3)

Noting that no direct observations of these stresses are available, the main problem in this paper is estimation of the stresses. A further limitation on such estimation involves the degree of lubrication at the base of the ice. This feature is crucial to understanding glacial dynamics and predicting the glacier’s response to climatic controls. However, this feature is not observable directly.

Treating the flow parameter A as a constant, the surface velocity u(x) is approximated by

(4)

where u _bx is sliding velocity and n is a flow parameter (e.g. Reference PatersonPaterson, 1994, p.251, equation 21). We fix n = 3 in this paper. Motivated by Equation (3), we model τ _bx as

(5)

where η_x is a corrector process needed to explain the velocity data. To clarify, we substitute Equation (5) into Equation (4) and estimate all quantities in light of the data. If η_x is found to be negative in some region, the velocities in that region are smaller than those expected based on a simple balance between driving and basal shear stresses. Hence, such η_x indicates additional resistive stresses acting on the flow. Alternatively, positive η_x suggests that the glacier is moving faster than expected, perhaps indicating an area of high lubrication at the base.

The simple interpretations above are of course flawed to some degree; we are attempting to use a single quantity η_x to explain the combined effect of at least two phenomena. For example, there may well be regions of both high lubrication and high resistive stress due to lateral effects. These are confounded in the simple approach here. This also suggests a confounding in modeling both sliding velocity u _bx and η_x as relatively unconstrained functions of x. We can break this confounding either by introducing very strong priors to constrain these processes or by constraining one of the processes while leaving the other relatively free and responsive to the observations.

In this paper, we compromise between these two notions. First, no matter how the modeling is carried out, we found that the data mandated the use of at least two different models. The estimated value of this change point between models is indicated in Figure 1. This estimate is close to an observable lineament. As will be indicated, creep flow appears to be the dominant mechanism upstream of this change point. Near the change point, sliding gains critical importance. We note that Reference Joughin, Fahnestock, MacAyeal, Bamber and GogineniJoughin and others (2001) suggest that there is a till plain downstream of where our analysis ends. We hope our results cast light on the nature of the bed between our change point and the start of the till plain.

Fig. 1.

Fifty posterior realizations of velocity profiles and the posterior mean (black) of velocities based on 2000 realizations. The original velocity data are shown (black dashes). Ensemble members are color-coded by importance-sampling Monte Carlo probability weights α (grey: α > 0.01; light grey: α < 0.01).

Based on the previous notions, we specified sliding velocity as a constant upstream of the change point, and used a relatively uninformative prior for η_x (see section 3.4). Downstream of the change point, we used approximations to Weertman-type models for sliding velocity (e.g. Reference PatersonPaterson, 1994, ch. 7) of the form

(6)

for various combinations of p and q with τ _bx given in Equation (5), i.e. without a corrector process. We found that the data strongly favored the selection of p = q = 1; we assume these values throughout the rest of this paper. Note that with these specifications, Equation (6) reduces to

(7)

We also tried using Equation (4) with no corrector process, no change point and u _bx as given in Equation (6) for all x and various combinations of p and q. None of these models was competitive with the model presented here.

Another interesting issue arises in the inference of driving stress based on observations of s and H. Reliance on the slope of the upper ice surface in Equation (2) implies that results are very sensitive to small-scale variations in surface topography and to small-scale, perhaps unimportant, variations in ice thickness. Hence, driving stress is usually estimated based on averaging over horizontal distances of a few ice thicknesses to eliminate small-scale flow features not important to the flow (Reference Kamb and EchelmeyerKamb and Echelmeyer, 1986; Reference PatersonPaterson, 1994). Indeed, if averaging is not done, driving-stress estimates exhibit unreasonably large variations. Hence, we incorporate a statistical smoothing step in our Bayesian model.

1.3. Data

We assembled surface topography and ice-thickness observations for a portion of the NEGIS (Figs 2 and 3). The data were gathered as part of the Program for Arctic Regional Climate Assessment (PARCA). Surface topography and ice thickness were sampled every few hundred meters using equipment mounted on the Wallops Flight Facility P-3 aircraft. Surface velocity data were calculated (personal communication from I. Joughin, 2002) dataset. The three derived datasets are S (surface topography), B (basal topography), and U (surface velocities). Datasets S and B are shown in Figure 2; U in Figure 1.

Fig. 2.

Observations of surface elevation and basal elevation along a 400 km profile of the NEGIS.

Fig. 3.

The basal elevation (black dashes), 50 realizations (grey) of smoothed basal topography from the posterior distribution conditional on those data, and the corresponding posterior mean (black) based on 2000 ensemble members.

1.4. Hierarchical Bayesian analysis

One rationale for the Bayesian approach is that it is a mathematically rigorous method for combining information in the presence of uncertainty. Two primary sources of information are available for inferring the unknown quantities of interest: (1) observations or data that convey some information regarding those unknowns, and (2) prior information based on scientific reasoning regarding the unknowns, including physical models as well as past experience and data. Both sources of information are subject to uncertainties. In the Bayesian approach, all such uncertainties are modeled probabilistically. The primary computational tool used is Bayes’ theorem, which is simply a result in probability theory relating various conditional distributions. However, the Bayesian view of statistical modeling and analysis involves more than the simple application of probability theory. It is a paradigm that involves the modeling of unknowns as random variables and using observations to update that modeling effort.

Many procedures used in the analysis of geophysical data are Bayesian or approximately Bayesian. For example, the Kalman filter is a Bayesian procedure (e.g. Reference Meinhold and SingpurwallaMeinhold and Singpurwalla, 1983). Other examples associated with data assimilation are discussed by Reference LorencLorenc (1986) and Reference Evensen and van LeeuwenEvensen and Van Leeuwen (2000). A general review of Bayesian analysis intended for geophysical audiences is given by Reference Wikle and BerlinerWikle and Berliner (2006); see also Reference EpsteinEpstein (1985) and Reference TarantolaTarantola (1987).

Nevertheless, the full power of the methodology has only recently been making progress in geophysics (e.g. Reference Berliner, Wikle and CressieBerliner and others, 2000; Reference Wikle, Berliner and MilliffWikle and others, 2003; Reference Tebaldi, Smith, Nychka and MearnsTebaldi and others, 2005). Advances in computation, through variants of Markov chain Monte Carlo (MCMC) algorithms, now enable hierarchical Bayesian modeling that is capable of dealing with the complexities in models and data that arise in geophysics. A skeleton of hierarchical reasoning begins with consideration of three basic collections of variables to be modeled: data (our observations) labeled y; process (those physical, state variables of interest (e.g. velocities, stresses, etc.)) labeled x; and parameters labeled θ, including unknown physical constants and parameters introduced in the statistical components of the model. We use the following notation for probability distributions: (1) the distribution of a random variable, say X, is written as [x]; and (2) the conditional distribution of Y given X = x is written as [y|x]. A review of the basic strategy and associated computations is given by Reference Gelman, Carlin, Stern and RubinGelman and others (2004).

Hierarchical thinking suggests a model with three primary components (e.g. Reference Berliner, Hanson and SilverBerliner, 1996):

1. 1. Data model: [y|x, θ];

2. 2. Prior process model: [x|θ]; and

3. 3. Prior on parameters: [θ].

Although the construction of these components is often not easy, once the modeling process is complete, Bayes’theorem yields the posterior distribution [x, θ|y] given by

(8)

The denominator [y] is the marginal distribution of the data. It can be viewed simply as that constant which normalizes the numerator ensuring that the posterior is indeed a probability distribution, depending upon the observed data but not the unknowns.

In the Bayesian treatment of parameters, even unknown constants are treated as if they are random. This is important since the definitions of some model parameters may be based on approximations. For example, flow parameters are not really unknown physical constants but rather functions of other unmodeled process variables (e.g. temperature or pressure) which introduce uncertainty. Hence, while the Bayesian approach may appear complicated compared to others, its advantage is that it quantifies and manages uncertainties through to final conclusions.

2.1. Translating physics into statistical models

In response to the issues introduced in section 1.2 regarding averaging of surface and thickness data for the estimation of driving stress, our first step is to model the basal topography b and the surface s as smoothed versions plus local noise, denoted generically by N. We therefore assume

(9)

where ψ_b (x) and ψ_s (x) are smooth versions of the base and surface, respectively.

Applying these definitions in Equation (2), we model smoothed driving stress as

(10)

where the smoothed ice thickness is . In all computations, we set ρ = 911 kg m⁻³ and g = 9.81 m s ^{− 2}.

We next define the partially smoothed basal shear stress in the region upstream of the currently undetermined change point, by applying Equation (5) for velocities taken in the negative x direction as

(11)

We say is only partially smoothed because we do not smooth the corrector process. As clarified later, η(x) is modeled as an unknown stochastic process. This enables further interpretation of the posterior behavior of η(x) as an indication of regions where successful or unsuccessful smoothing has occurred.

We then incorporate the velocity model Equation (4) (recall that we set n = 3) and the approximate sliding velocity model Equation (7) into a statistical model:

(12)

(13)

where k, A ₁, u_b , A ₂ and c are treated as unknown, random parameters and N_u (x) is an error process.

Before proceeding with the specifics, we emphasize that these formulae motivate probability models that are updated in light of the observations. More directly, we do not produce single estimates of smoothed processes ψ_s (x) and ψ_b (x) based on the data and substitute those estimates into the stress and velocity models. Rather, we produce an ensemble of these quantities simulated from their posterior distributions and propagate them forward into our models. This produces an ensemble of stresses and velocities, which in turn are updated based on the velocity observations.

2.2. Underlying probability theory

Our main task is the development of the following probability distributions:

where θ denotes the collection of all model parameters. Our goal is to obtain the posterior distribution [ψ_b , ψ_s , u, η, θ |B, S, U], which can then be used to obtain the posterior distribution of stresses.

Our main assumption regarding the data model is that the three primary datasets are conditionally independent given the true processes they represent and given the model parameters. Notationally, the data model takes the form

(14)

where notation such as θ _B is used to indicate those parameters (subsets of θ ) explicitly appearing in the indicated models. A possible objection to Equation (14) arises since the basal data B are actually computed as the difference between surface and thickness observations. The assumed conditional independence therefore may not hold. We checked our posterior results for indications of degrees of departure from this assumption and found none that would affect our results.

Our modeling of the processes begins with a probabilistic equality (i.e. this is not an assumption, but a fact):

(15)

Similarly, we have

(16)

We make the following two assumptions.

1. 1. In formulating the second term on the right side of Equation (15), we assume that the base b and the surface s are independent, conditional on their smoothed versions and the model parameters. We then obtain

   (17)

   

   where again notation such as θ _b is used to indicate appropriate subsets of θ . It is critical to note here that we are not assuming that the base and the surface are independent. Our modeling of both the base and surface is conditional upon smooth processes ψ_b and ψ_s . Our assumption is that the small-scale departures from those large-scale processes are independent.

2. 2. Consider the conditional model for u, namely the first term on the right side of Equation (16). We assume that the velocity profile depends on the base and surface only through their respective smoothed versions ψ_b and ψ_s , i.e.

   (18)

   

   The specifications of prior distributions for the parameters are given in Appendix A.

3.1. Basal model

First, the process model is described. We chose to use wavelets to model the smoothed basal topography ψ_b (x) because of their flexibility in representing highly variable processes and because we can easily control their smoothness.

Wavelets work best for equally spaced data where the number of data points is an integer power of 2. Hence, we partition the domain of the data into 2¹¹ = 2048 bins of equal length (189.5 m). Let denote the 2048-dimensional vector constructed by averaging b within each bin. Note that is not observed.

We used multi-resolution Daubechies wavelets (e.g. Bruce and Gao, Reference Bruce and Gao1996; Vidakovic, Reference Vidakovic1999). Two collections of wavelets were used. The first set is thought of as the smooth signal, while the second represents the detail signal. We consider a smooth signal of four wavelets. To this component we added a detail signal containing 28 more wavelets. All 32 wavelets are mutually orthogonal.

After converting to a discrete wavelet form, we obtain a linear model for :

(19)

where W is the 2048 × 32 matrix of discretized wavelet basis functions, C is the 32-dimensional vector of wavelet coefficients and Σ(ϕ ₁, ϕ ₂) is the correlation matrix of an autoregressive process of order 2 (AR(2)) with variance . Notationally, ψ_b = WC and . The selection of an AR(2) process to account for spatial dependence among the model errors (i.e. local variations in basal topography) was based on preliminary data analysis and practicality; our Bayesian computations require repeated inversion of a 2048 × 2048 matrix involving the inverse of Σ(ϕ ₁, ϕ ₂), which is straightforward for an AR(2) process.

We performed analyses for four resolutions, corresponding to 8, 16, 32 and 64 coefficients in Equation (19). We found that the model with k = 32 coefficients provided the best results in terms of maintaining fidelity to the thickness data (i.e. not over-smoothing), yet providing good inferences for velocity and therefore stress.

Turning to the data model, we define the basal data vector of length 2048 with the ith element given by the simple arithmetic average of those basal observations lying in bin i. Let n_i (i = 1, …, 2048) be the number of data points lying in bin i. In our case, none of the n_i were equal to zero; most were either 1 or 2.

Our data model (i.e. [B |b, θ _B ] in Equation (14)) is

(20)

where represents measurement-error variance and is a 2048 × 2048 matrix with diagonal elements equal to and off-diagonal elements equal to 0 (i.e. the elements of are assumed to be conditionally independent). This suggestion arises from a familiar result from elementary statistics. The sample mean of n independent observations with common mean and common variance σ ² has variance σ ² /n.

Specifications of priors on parameters are given in Appendix A.

3.2. Surface model

Our modeling strategy for [s|ψ_s , θ _s ] in Equation (17) separates the large-scale and small-scale behaviors of the surface. Recalling Equation (9), we suppose

(21)

where the large-scale surface is given by a function ψ_s , assumed known up to a low-dimensional set of parameters, and N_s is a zero-mean spatial stochastic process described in Appendix A. To model ψ_s , we rely on an analysis given in Paterson (Reference Paterson1994, p.243, equation 8) and assume the basic model for the surface:

(22)

where θ _s = (μ, K, L) are treated as unknown parameters. Although we set n = 3 here, we could also model n as an unknown.

We use only the large-scale surface Equation (22) to compute ice thickness, the surface derivative and hence the stress in Equation (10). Nevertheless, the presence of N_s is important in determining the data model in Equation (14). Under the strategy which uses Equation (22) to obtain the stress, we need a data model [S |ψ_s , θ _s , θ _S ], where θ _S includes measurement-error variances and parameters of the distribution of N_s in Equation (21); our specification is given by Equation (A3) in Appendix A.

3.3. Velocity model

For the data model [U |u, θ _U ], we assume that, conditional on the true velocities, the elements of U have independent Gaussian distributions with means equal to the true velocities at the corresponding locations and common variances . That is, at each observation location x, we assume that

(23)

where m(x) is a zero-mean, Gaussian measurement error with variance . Substituting for u(x) as prescribed in Equations (12) and (13), we have

(24)

(25)

We next combine the model errors N_u and traditional measurement errors m into a single error

(26)

assuming that e_U (x) are independent, zero-mean, Gaussian random variables with common variance . (Note that these steps mean that there is no need for further use of the symbol u.) This simplification can be relaxed, but would require additional modeling of the process N_u .

To simplify the presentation of the statistical model, we write our models in vector notation. Recalling Equations (9–11) and (19), we define the vector of smoothed ice thicknesses

(27)

where ψ _s is the 2048-dimensional vector of smoothed surface elevations. Similarly, we define the vector of smoothed values of driving stress τ _d as

(28)

where * indicates the Hadamard product, i.e. we compute the vector of element-wise products of the coordinates of the smoothed thickness and the derivative of the smoothed surface.

Next, we let η denote the vector of stress corrections at the locations of velocity observation and set

(29)

From Equation (12) we model u, the vector of true velocities at the observation locations, as a linear function of the corresponding coordinates of smoothed thickness multiplied by the third power of coordinates of τ _b.

In preliminary data analyses, we noted that at least two models (one for small x and another for large x) are needed. Let x = c be an unknown change point and consider different linear functions above and below the change point. Finally, the model for the velocity data vector U is

(30)

where subscripts 1 and 2 indicate the varying dimensions of the vectors (a vector with all elements equal to 1) depending on the value of the change point c and e _U which are measurement and model errors, respectively (recall Equation (26)).

3.4. Stress corrector process

Since the corrector process arises from unmodeled and unobserved effects (e.g. longitudinal and lateral drag, basal lubrication, local under- or over-smoothing, etc.), we assume a relatively uninformed prior for η . Specifically, our prior [ η ] indicates that each element has a uniform distribution on the interval (−100 kPa, 100 kPa) and that all elements are mutually independent.

3.5. Bayesian calculations

Although we can write down Bayes’ theorem for the posterior distribution of all unknowns conditional on the observations, the result is typically not computable in closed form. We use a Monte Carlo approach that produces an ensemble of realizations from the target posterior distribution.

The method relies on the emerging technology of MCMC. The idea of MCMC is to simulate a Markov chain that has been carefully designed so that its stationary distribution coincides with the target posterior distribution. It follows that, after a burn-in or transience period, the generated realizations of the chain comprise a simulated sample from the posterior. Data analysis (often known as ‘output analysis’) is performed on this sample to produce the desired inferences.

In our case, direct use of MCMC is possible but challenging, primarily due to the non-linearities present in Equations (3) and (12). Hence, we combine MCMC with the technique of importance-sampling Monte Carlo (ISMC). Consider a setting in which direct simulation from a target distribution is difficult or inefficient. In ISMC, one generates an ensemble from another more manageable distribution. The theory of ISMC provides formulae for the calculation of weights that are used to re-weight the ensemble, permitting inferences relative to the original target. General introductions to both MCMC and ISMC can be found in Robert and Casella (Reference Robert and Casella1999). An illustration of these technologies in a geophysical problem is given in Berliner and others (Reference Berliner, Milliff and Wikle2003).

An outline of the calculations used here is as follows. We first run separate, independent MCMC algorithms (simple versions known as Gibbs samplers) for the basal model and the surface model. These runs produce ensembles from the posterior distributions and [ψ_s , θ _s , θ _S| S]. They are then used in conjunction with the velocity model [u|η, ψ_b , ψ_s , θ _b , θ _s , θ _u ] (recall Equation (18)) to simulate velocities conditional on and S. (As described in Appendix A, we actually use a portion of S.)

To incorporate the velocity data U, we re-weight all of these samples using ISMC results, yielding the ensemble used to summarize the full posterior distribution. Technical details regarding the actual implementation are presented in Appendix B.