Variational coupling of Plinian column models and data: Application to El Chichón volcano

Journal of Volcanology and Geothermal Research 175 (2008) 501–508

Contents lists available at ScienceDirect

Journal of Volcanology and Geothermal Research j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j vo l g e o r e s

Variational coupling of Plinian column models and data: Application to El Chichón volcano Isabelle Charpentier Laboratoire de Physique et Mécanique des Matériaux, UMR CNRS 7554, Ile du Saulcy, 57045 Metz Cedex 1, France

A R T I C L E

I N F O

Article history: Accepted 12 February 2008 Available online 16 April 2008 Keywords: modelling inverse problem data assimilation numerical simulation

A B S T R A C T Many efforts have been done in the modelling of Plinian columns. However, until now, the inverse problem of the reconstruction of a Plinian event from observed data has been only roughly tackled. This paper discusses the efficient variational data assimilation (VDA) that manages the optimization iterate sequence by means of gradient computations. Theoretical developments of VDA are presented for both the Woods equations and a generic system allowing for the modelling of the Navier–Stokes equations. Two- and three-dimensional unsteady Plinian models being based on the latters, VDA could be clearly performed on their differentiable subset of equations. Although no column data are available for the 1982 El Chichón eruptions, identical twin experiments were used to check the abilities of the VDA approach. VDA experiments coupled with a tephra tephra-fall model could be used to obtain isopach and isopleth maps of past eruptions. Performing VDA experiments with a Plinian model coupled to a tephra-fall model could be an interesting solution to exploit isopach and isopleth maps of the past eruptions of El Chichón. Some of the numerical experiments presented in this paper are downloadable in the PliniAD package which is freely available from the author's webpage. © 2008 Elsevier B.V. All rights reserved.

1. Introduction The formation of Plinian columns is a subject of fundamental importance in volcanology and climatology as these explosive eruptions are the mechanism by which large amounts of particles and aerosols are released in the high atmosphere. As discussed in Matson (1984), 1982 El Chichón major eruptions on March 29 (0532 UT), April 4 (0135 UT), and April 4 (1122 UT), as well as Mount St. Helens (18 May, 1980, USA) and Alaid Volcano (April 27, 1981 and May 1, 1981, Russia) eruptions, were the largest documented by operational weather satellite data since the beginning of such data collection in 1966. Nevertheless, evidences of a strong stratospheric penetration were noticed by GOES and NOAA-6 satellites for the El Chichón eruptions only (Matson, 1984). In that sense, modern studies on climate effects of Plinian eruptions have really begun with the 1982 El Chichón eruption and ash cloud satellite data (see Robock, 2000 for a review on volcanic eruptions and climate). They become more organized, involving much more satellite observations, with the November 7, 1991 Pinatubo eruption (Philippines). Plinian columns are hot mixtures of steam, volatiles, liquid water, and solid fragments ejected through the atmosphere at high velocities. In the thrust region, the motion of the column results from the high pressure difference (Nairn, 1976; Kieffer, 1981; Kieffer

E-mail address: [email protected]. 0377-0273/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jvolgeores.2008.02.021

and Sturtevant, 1981) between the mixture in the volcanic conduit and the atmosphere at the vent. The mixture, forced into the atmosphere, looses kinetic energy and ingests ambient air. Through entrainment and heating of air, this turbulent flow may become buoyant and reach altitudes that exceed the tropopause. Complex processes ranging from the microscales to the macroscales govern the column dynamics. First modellings were proposed three decades ago assuming a fine-grained pyroclast distribution for which the transfer of heat and momentum between particles and gas is rapid. Thus Wilson (1976), Wilson et al. (1978) and Settle (1978) based their approach on the mass-, momentum-, and energy-conservation equations of the motion of one fluid into another following the fluid mechanics theories of jets (Prandtl, 1949) and buoyant plumes (Morton et al., 1956). Sparks (1986) improved the earlier models by taking into account the influence of the magma temperature and the vertical atmospheric profiles of pressure and temperature. In a more sophisticated model, Woods (1988) joined together the jet and plume models to propose a “top hat formulation” that takes into account the conservation of thermal energy and the entrainment assumption introduced by Morton et al. (1956). A review of models dedicated to the general problem of plume formation and dynamics is presented in Sparks et al. (1997). Several 2-dimensional and 3dimensional unsteady models were proposed (Valentine and Wolhetz, 1989; Dobran et al., 1993; Oberhuber et al., 1998). Among them, ATHAM (Active Tracer High resolution Atmospheric Model) developed in Herzog et al. (1998) and Oberhuber et al. (1998) solves

502

I. Charpentier / Journal of Volcanology and Geothermal Research 175 (2008) 501–508

the unsteady non-hydrostatic Navier–Stokes equations. In this model, transported tracers influence the dynamics of the flow. Recently, scientific efforts were devoted to the design of complex models that deal with volcanic features ranging from microscale to macroscale phenomena. The joint work proposed in Textor et al. (2005) is probably impressive as it discusses models from the conduit processes to the meteorological and climatological effects of ash and aerosol dispersion. The reader is referred to Textor et al. (2005) for an exhaustive bibliography on the modelling of Plinian column and related processes. On the one hand, volcanological models attempt to capture the main aspects of complex processes in a mathematical formulation whose validity lies on its ability to yield results in accordance with observation. On the other hand both the models and the observation systems play a role in the evaluation and the prediction of environmental risks. Thus coupling models and data has become an important scientific goal in most of the geophysical communities. In these regards, the first question that emerges is about the sensitivity of the different state variables of a model to a given perturbation in the inputs. This is a keypoint since it reveals the stability of the system and the relative relevance of the different variables of the model. This question can be answered at different levels considering either a qualitative observation of results stemming from successive trials (as usually performed in volcanology), or quantitative methods based on the computation of derivatives as previously proposed in Charpentier and Espíndola (2005a,b) and Ishimine (2006). More generally, optimal control theory (Lions, 1971) and gradient computation have become important tools in scientific computing. Applications are sensitivity analysis, identification (Chavent, 1975), forecasting (Le Dimet and Talagrand, 1986; Courtier and Talagrand, 1987; Talagrand and Courtier, 1987; Lewis and Derber, 1985), and shape optimization (Pironneau, 1983). It should be stressed that an impressive number of articles derived from these pioneering works. To our best knowledge, the variational data assimilation (VDA) technique, introduced by Le Dimet and Talagrand (1986) and Lewis and Derber (1985) in the context of meteorology, has never been used within the Plinian column context before. VDA is devoted to the solution of a minimization problem combining quantitative information on a phenomenon with modelling results by means of the optimal control theory. Such an inverse problem makes use of an objective function that measures some discrepancies between acquired data and numerical solutions. A minimization problem involving a compound function (objective function + model) is set up for the identification of some inputs of the model (hereafter named controls) from which the solution fits at best the observed data. Assuming that both the model and the objective function are differentiable, VDA proposes to solve this inverse problem looking for a control that cancels out the gradient of this compound function. Great advances have been achieved with the development of automatic differentiation (AD) tools (Griewank, 2000) that may be applied to large numerical codes (Charpentier, 2000). Nowadays VDA has successfully succeeded to statistical approaches for weather forecast, being daily used in most of the operational meteorological centers. Since then, VDA has been widely applied to geophysical flow models in oceanology, hydrology, and atmospheric chemistry. Many other inverse problems may be tackled with VDA as soon as data are available and relevant. The first historical report of a past volcanic activity of El Chichón volcano was published by Canul and Rocha (1981) a few months before the 1982 eruption. During the El Chichón Plinian eruptions (1982), data were recorded by the GOES-East weather satellite in the visible and IR bands (National Earth Satellite Service (Ness), 1982; Matson, 1984) whereas meteorological stations at Veracruz and Merida (380 km and 590 km northwest and northeast of El Chichón respectively) measured wind profiles by means of radiosondes. The major eruptions produced widespread tephra-fall deposits and

generated dramatic pyroclastic surges and flows (Sigurdsson et al., 1984). In their paper, Carey and Sigurdsson (1986) presented isopach and isopleth maps and used a tephra-fall model (Carey and Sparks, 1986) to perform inverse experiments. More recently, Espindola et al. (2000) presented a revised composite stratigraphic column, showing that El Chichón erupted at least 12 times during the last 8000 years including the 1982 eruption. Of these eruptions only the 550 year B.P. eruption produced widespread tephra-fall deposit. Unfortunately, as mentioned by Sparks et al. (1997), a few Plinian eruptions have been filmed with the purpose of making scientific measurements. In particular, no scientific ground-based images of the column were filmed for 1982 El Chichón eruptions (that did not occur at daylight). Thus, in this paper, the interest of VDA is brought to the fore considering identical twin experiments (Charpentier, 2007) that involve synthetic data generated from the model (Bunge et al., 2003). White Gaussian noises are then added to mimic actual data. Such a study is attractive since modelling and observation errors are under “control”. Moreover, as any variable of the model is numerically observable, it allows for an exhaustive investigation of the optimization capabilities. Twin experiments are of primary importance: if the inversion process is not convenient in this numerical context, it will generally not be more convenient using actual data. First numerical experiments based on twin experiments were discussed in Charpentier (2007) to prove the interest of ground-based observations of the tropospheric part of the Plinian columns. In the same paper, VDA was also applied with success for the identification of the entrainment function as cubic splines considering actual data measured by Calder et al. (1997) from images of growth profiles (Hoblitt, 1986). The PliniAD package discussed in the paper is freely downloadable from the author's webpage1. It proposes a perfect twin experiment, as well as the original model and two automatically differentiated codes, written within the El Chichón volcanic context. In the paper, VDA is applied to the simultaneous identification of the inflow conditions and the entrainment parameter performing several twin experiments. The layout of the paper is the following. In Section 2 Woods (1988) equations are modified to be in the VDA general framework introduced by Le Dimet and Talagrand (1986). Section 3 discusses the differentiation of the model from both theoretical and numerical points of view. Numerical VDA results are presented in Section 4. Some prospective elements are proposed as a conclusion. 2. One-phase 1-dimensional modelling In this section, we shall use the model developed by Woods (1988) since it is the most general among the single-phase models. Once converted into an Ordinary Differential Equation (ODE) system, this model may be written under the generic form introduced by Le Dimet and Talagrand (1986) in their presentation of VDA for Partial Differential Equation (PDE) problems. Originally this generic model was designed to the inverse modelling of the meteorological PDE flows, but the VDA framework is convenient for other geophysical flows, including 2- and 3-dimensional PDE models of Plinian columns. 2.1. Original equations The original equations of Woods (1988) are written in terms of velocity u, radius r, bulk density β, closely related to gas fraction n, and temperature θ of the column. They make use of physical parameters such as volcanic gas constant Rg0, bulk specific heat Cp0 of mixture in column, and density of pyroclasts σ, mean variables Rg and Cp being estimated along the column. Woods (1988) equations operate from inflow conditions u0, r0, n0 and θ0 at the vent (z = z0). 1

lpmm.fr/charpentier/pliniad.

I. Charpentier / Journal of Volcanology and Geothermal Research 175 (2008) 501–508

503

Using these notations and omitting the dependence on height z, equations are: (1) mass balance along the column (z0, ztop):

Related boundary conditions are deduced through the bijective change of variable T that maps (u0, θ0, r0, n0) to (u0, E0, F0, n0):

dbur2 ¼ 2uera; dz

T ðu0 ; h0 ; r0 ; n0 Þ ¼




ð1Þ

where α is the density of the ambient air, ε is a parameterized entrainment function representing the average form in which air is trapped in the column, and ztop is the height of the column, (2) momentum balance along the column: du ða  bÞr 2 u ¼g  ð2ueraÞ; dz bur2 bur2

ð2Þ

where g is the acceleration of gravity, (3) energy balance along the column:

 dCp h 2uera u2 g ða  bÞr 2  uð2ueraÞ ¼ gu Ca T  Cp h  ; 2 2 dz 2 bur bur

ð3Þ

where Ca is the specific heat of air and T the ambient air temperature. These equations are completed by functional relationships between variables and parameters:

 8   1n 1 1 nRg h n0 > > ; Rg ¼ Ra þ Rg0  Ra ; < ¼ ð1  nÞ þ P b r n 1  n0 2   ð1  nÞ b u r > > : n ¼ 1 þ ðn0  1Þ 0 020 ; Cp ¼ Ca þ Cp0  Ca ; ð1  n0 Þ bur

ð4Þ

where Ra is the gas constant of the air. The pressure P may be derived from temperature T from standard atmospheric models. 2.2. ODE system The former system of equations is turned into an ODE system of 4 equations involving 4 unknowns by introducing the mass flux variable F = βur2 and the energy variable E = Cpθ. Two additional variables D = (α− β)/(βu) and G = 2εαur/(βur2) are used to simplify the equations. The balance equations may be written as: 8 dF > > ¼ FG; > > dz > > > < du ¼ gD  Gu; dz > >   dE > > ¼ G Ca T  E  u2 =2  g  uð gD  GuÞ > > > : dz ¼ GðCa T  EÞ  g  gDu þ Gu2 =2:

 ð8Þ

In the sequel, the notations we use are mainly those proposed in Ide et al. (1997). The Plinian ODE system is summarized as: ( ðMÞ

  du u ¼ Mðu Þ in z0 ; ztop ; dz u ðz0 Þ ¼ u 0 ;

ð9Þ

where M is the nonlinear operator that describes the dynamics of the model, and u0 = (u0, E0, F0, n0) is a boundary condition vector belonging to the set U0 of admissible boundary conditions for a Plinian column. For VDA purposes, operator M is assumed to be differentiable with respect to the boundary conditions. Such a generic system allows for the modelling of spatio-temporal PDE problems such as the Navier–Stokes equations. As 2-dimensional or 3-dimensional unsteady Plinian models are based on the latters, the theoretical developments presented below may be applied to the differentiable part of their equations. The Plinian code related to Eqs. (5)–(8) is part of the PliniAD package. The latter may be freely downloaded from the author's webpage1. (see footnote 1). The Plinian code is made of 5 routines and an include file containing data (volcanic, atmospheric, …) and common variables for the observed data management. The main routine, Plinian, initializes inflow conditions u0, theta0, r0, n0 to user-defined values and the objective function result j to 0 (see Section 2.3 for details) before a call to the model calM. State variable values along the column are stored in a file for drawing and observation purposes. Following the LaTeX calligraphic letter command, the calM routines organize the solution of M through a handcoded RK4 scheme. During each step of the latter, 4 calls to routine M allowing for the evaluation of both the M operator present in the right-hand side of Eqs. (5)–(6) and the functional relationships of Eq. (7). The objective function result calJ, is updated at the end of each RK4 step. Standard atmosphere variables are computed in M (and once in calM). The call graph is the following one.

ð5Þ

Differentiating the gas fraction Eq. (4) with respect to z allows for the construction of the last ODE: dn d dF=dz F0 G ¼ ð1 þ ðn0  1ÞF0 =F Þ ¼ ðn0  1ÞF0 : ¼ ð1  n0 Þ dz dz F F2

u0 r02 ; n0 ð1  n0 Þ=r þ n0 Rg0 h0 =P0 ¼ ðu0 ; E0 ; F0 ; n0 Þ: u0 ; Cp0 h0 ;

ð6Þ

Obviously variables θ, r, β, D, G and parameters Rg and Cp may be deduced from the main variable vector (u, E, F, n) from the following ordered sequence of equations: 8
   n0 1  n   ð1  nÞ > > R ; ¼ R þ R  R ; Cp ¼ Ca þ Cp0  Ca > g a a g0 > > 1  n n ð1  n0 Þ 0 > > < h ¼ E=Cp ;
 1 nRg h 1 > > > b ¼ ð1  nÞ þ ; > > P > pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r >   : r ¼ F=ðbuÞ; D ¼ ða  bÞ=ðbuÞ; G ¼ 2eaur= bur2 :

2.3. Inverse modelling The aim of any data assimilation process is the identification of some inputs of a model, also named control, by including observed data vo (belonging to the set U o of acquired data) in the modelling. Of course, this model is assumed to be able to compute realistic states of the system under observation when provided with adequate inputs (initial or boundary conditions, parameters…). The noises resulting from observation are assumed to be white Gaussian ones. Data are introduced by means of an objective function J :U → IR J ðuÞ ¼

ð7Þ

1 2

Z

ztop

z0

jj Hu  vo jj 2 ;

ð10Þ

where H is the observation operator that maps the set U of admissible states of the column to the set U o of acquired data.

504

I. Charpentier / Journal of Volcanology and Geothermal Research 175 (2008) 501–508

For the sake of clarity, let us first discuss about the identification of vector u0 in the generic framework. One searches a solution um 0 of Eq. (9) that fits at best data vo along the column. That is: Find u m 0 a U0

  such that J &M u m u0 aU: 0 VJ &Mðu 0 Þ 8u

ð11Þ

Assuming the compound function J &M is locally convex and differentiable with respect to the control u0, the VDA approach allows for the solution of the minimization problem (11) as: Find u m 0a U

such that jðJ &MÞðu 0 Þ ¼ 0:

ð12Þ

One notices that the former two hypotheses are satisfied as soon as the nonlinear operator M of system (9) is a differentiable one. In the Plinian case, we are interested in the identification of the inflow conditions (u0, θ0, r0, n0). The compound function to minimize is actually J M T . VDA may be used since the change of variables (8), the (Woods, 1988) Eqs. (5)–(7) and the objective function (10) are clearly differentiable. Explanations about the numerical computation of the gradient ∇(J M)(u0) are provided in the next section, pointing out automatic differentiation software abilities.

& &

&

3. Gradient computation Gradient computation becomes an important tool in scientific computing. Applications are sensitivity analysis (Charpentier and Espíndola, 2005a,b; Ishimine, 2006), parameter identification (Chavent, 1975), reconstruction of events and forecasting (Le Dimet and Talagrand, 1986; Courtier and Talagrand, 1987; Talagrand and Courtier, 1987; Lewis and Derber, 1985), and shape optimization (Pironneau, 1983). The two methods for getting differentiated codes are: writing the code from the derivatives of the continuous mathematical equations, or differentiating the code that discretizes the continuous mathematical equations. Then, for each method, two modes of differentiation may be considered. The following two subsections describe the tangent linear and adjoint methods for the computation of the gradient of J appearing in Eq. (12). Automatic differentiation is then briefly presented as it avoids the tedious task of a handwritten derivation of the codes. 3.1. Forward mode of differentiation The tangent linear differentiation of J ∘M is based on a mechanical use of the chain rule, that is: ju 0 ðJ &MÞðu 0 Þ  du u0 ¼ ju uJ ðMðu 0 ÞÞ  ju0 Mðu 0 Þ  du u0 ;

ð13Þ

where the Jacobian ∇uJ (M(u0)) is the matrix of derivatives of J calculated with respect to variable u and evaluated at point M(u0). Variable δu0 ∈ U 0 is a given direction of perturbation (of the boundary condition u0). Applied to the generic problem (9)–(10), one deduces that:

(

    ddu u AM ðu Þ  du ¼ u in z0 ; ztop ; ðMl Þ dz Au u u0 ; du uðz0 Þ ¼ du Z ztop u0 ¼ H⁎  ðHu u  vo Þ  dudz; jðJ &MÞðu 0 Þ  du

ð14Þ

z0

where the tangent linear perturbation variable δu ∈ U is computed along the column through the solution of M1. The tangent linear model is evaluated along the trajectory indicated by the presence of variable u in the right-hand side terms of Eq. (14). From the computational point of view, the trajectory is provided through the solution of the generic problem (9)–(10). One observes that the computation of all the components of ∇J requires as many evaluations of Eq. (14) as the dimension p of the control u0. According to theoretical bounds computed

by Morgensten (1985), the CPU cost of the gradient computation is at most equal to 4p times the CPU cost of the original code integration. Although trivial rules are used, the differentiation of the continuous equations is usually a tedious task. This may be observed when differentiating the Plinian Eqs. (5)–(6) with respect to the 4 control variables u0, θ0, n0 and r0. The direction of perturbation is (δu0, δθ0, δn0, δr0). Following the chain rule, any variable presenting a dependence with respect to one of the controls is active for the differentiation. Thus the differentiation of the change of variable T implies the computation of: 8 < db ¼ b20 ðdn0 =r  Rg0 ðdn0 h0 þ n0 dh0 Þ=P0 Þ; db ¼ db0 u0 r02 þ b0 du0 r02 þ 2b0 u0 r0 dr0; : db ¼ Cp0 d0 :

ð15Þ

Parameters Cp0 and Rg0 are not differentiated because they are user-defined constants at vent. On the contrary, parameters (Cp and Rg) calculated along the column are active ones. The tangent linear Plinian ODE system is: 8 ddF > > ¼ dFG þ FdG > > > dz > > > ddu > < ¼ gdD  dGu  Gdu; dz ddE > > > ¼ dGðCa T  EÞ  GdE  gdDu  gDdu þ dGu2 =2 þ Gudu; > > dz > > > > : ddn ¼ ð1  n0 ÞF0 dFG þ FdG : dz F2

ð16Þ

where, following the chain rule, tangent variables δD and δG have to be known too. One has for instance: dD ¼ ½dbðbuÞ  ða  bÞðdbu þ bduÞ=ðbuÞ2 ;

ð17Þ

that itself requires the knowledge of several derivatives: 8     > dR ¼  Rg0  Ra n0 =ð1  n0 Þ dn=n2 ; > > < dC ¼ dnC  C =ð1  n Þ; p p0 0  a dh ¼ dECp  EdCp =Cp2 ; > >  

> : db ¼ b2 dn=r  dnRg h þ ndRgh þ nRg dh =P :

ð18Þ

The differentiation of G and r is not reported here for the sake of brevity. Moreover, such a computation, which is essentially of mathematical interest, is not really needed for numerical purposes since efficient AD tools allows for the differentiation of complex computer codes. 3.2. Reverse mode of differentiation The computation of the gradient components may be solved upside down. From a theoretical point of view this may be explained as in Talagrand and Courtier (1987). Let us introduce a scalar perturbation δω = ∇(J ∘M)(u0) δu0 of the discrepancy ω = J ∘M(u0) ∈ IR and 2 inner products b.,.N and [.,.] on IR and U 0 respectively. Up to a scaling, one may choose δω = ||δω|| = 1. One deduces: u0 N 1 ¼ bdx; dxN ¼ bdx; jðJ &MÞðu 0 Þ  du ¼ dx⁎  jðJ &MÞðu 0 Þ  du u0 ¼ hðjðJ &MÞðu 0 ÞÞ⁎ dx⁎  du u0i ˆ ðz0 Þ; du ¼ ðjðJ &MÞðu 0 ÞÞ⁎ dx; du u0 ¼ ½ u u0 ;

ð19Þ

&

where the adjoint variable û(z0) = (∇(J M)(u0))⁎ belongs to the space of controls U 0. In their original proof, Le Dimet and Talagrand (1986) used a Lagrangian approach to construct the adjoint model (20) from the generic problem (9)–(10): 8  ⁎ > > ˆ þ AM ðu Þ u   > :u ˆ ztop ¼ 0;

  in z0 ; ztop ;

ð20Þ

I. Charpentier / Journal of Volcanology and Geothermal Research 175 (2008) 501–508

ˆ∇ enables the computation of the gradient as: which solution u jJ •Mðu 0 Þ ¼  u ˆ ðz0 Þ:

ð21Þ

From a computational point of view, the adjoint code is integrated backward from the top of the column to the vent. This signifies that the trajectory of evaluation has to be cautiously taken into account in the implementation of any adjoint code. Following Eq. (21) the full gradient requires a unique evaluation of the adjoint model: the CPU cost is independent from the dimension of the control space. Theoretical results (Morgensten, 1985) indicate that, for any assignment statement, the ratio between the number of operations of the adjoint statement and the number of operations of the original one is lower than 4. As a result, the adjoint code enables to overcome the deadlock of the large dimension of u0 as it occurs in weather forecast experiments for instance. Note that the adjoint derivatives produced by the two methods may differ depending on the equations and their discretization. In such case, a theoretical analysis is necessary to choose between the two (Roche, 1997). In peculiar, such a difference is observed when differentiating in the adjoint mode a Leap-Frog scheme used for advancing in time the solution of a Navier–Stokes PDE system (Sirkes and Tziperman, 1997). In this case, the differentiation of the computer code has to be preferred as it allows for gradient computation ∇(J M) that are in accordance with the numerical computation of J M. Adjoint equations of the Plinian model are not provided in the paper since this tedious task may be performed using an AD tool.

•

•

505

The atmosphere routine is not differentiated since it does not depend on the independent (control) variables u0, θ0, r0, n0. The call graph of the adjoint code generated by Tapenade is too large to be fully displayed. Indeed, the trajectory management is performed by means of specific routines that store variable values during a preliminary computation of the Plinian model, and restore them during the derivative evaluation. The Plinian routines of the call graph of the adjoint model calM_B are:

Following Eq. (19), adjoint routines are called in the reverse order because:  u ˆ ðz0 Þ ¼ ðjðJ •MÞðu 0 ÞÞ⁎ ¼ ðjM ðu 0 ÞÞ⁎  jJ ðMðu 0 ÞÞ⁎ :

ð22Þ

3.3. Automatic differentiation The first differentiated codes were hand-coded using well known rules such as “the derivative of a sum is the sum of the derivatives”, the chain rule being applied to statements containing several operators. The differentiation becomes a fastidious and error prone task when large codes are considered. This is even truer in the adjoint mode of differentiation as it requires a complex management of trajectory. Nowadays, powerful automatic tools allow for the differentiation of general codes representing mathematical functions (Griewank, 2000). The reader is referred to Internet site autodiff.com for an exhaustive presentation of the existing softwares. The PliniAD package provides both the tangent linear and adjoint codes generated with the Tapenade on-line automatic differentiation engine (Hascoët and Pascual, 2004) from the Plinian model provided in the package. In practice, the differentiation step may be performed with Tapenade by: (1) uploading the sources files and the include file, (2) choosing the top routine name (calM), (3) choosing the independent variables (u0, theta0, r0, n0), (4) differentiating in Tangent Mode, or in Reverse Mode (adjoint mode), (5) downloading the generated files. In the automatically differentiated codes, tangent linear and adjoint variable names are respectively suffixed by d and b, whereas routine names are suffixed by _D and _B. The interested reader may compared the automatically generated tangent linear statements of the M_D subroutine with the tangent linear Eqs. (17)–(18). M_D also contains the original Plinian equations in order to provide a correct trajectory for the evaluation of the derived statements. The call graph of the tangent linear code is:

With both codes, derivatives are computable up to the machine precision, which is of prime interest for the solution of inverse problems requiring accurate gradient evaluations. The PliniAD package also provides a Taylor test for each of the two differentiated codes to verify the correctness of the differentiated codes. The correctness of the tangent linear code is checked through a comparison of Eq. (23) of its gradient computation with a first order Taylor approximation: i rw

¼

  jJ oM u0 þ xei  J oMðu0 ÞÞj xjjðJ oMÞ⁎ ðu0 Þ  ei j

; 8ei aIR4 :

ð23Þ

The test presented in Table 1 is correct since the ratio rω tends linearly towards 1. This is true up to ω≃10− 5 which is the optimal value for the finite difference approximation of temperature component of the gradient. Afterwards, the subtraction of too close floatingpoint numbers (very small steps ω) leads to a large cancellation error that dominates the truncation error in the approximation of the gradient by the finite difference scheme. Same results are obtained for a gradient computed with the adjoint code. 4. Numerical results 4.1. Sensitivity analysis Given inflow conditions and a direction of perturbation, the tangent linear code allows for straightforward sensitivity computation. It

Table 1 Result of the Taylor test for a perturbation δθ0 = ωK of the boundary temperature ω

rω

ω

rω

1.E + 01 1.E + 00 1.E − 01 1.E − 02 1.E − 03 1.E − 04

0.920678496 0.992001183 0.999199444 0.999919937 0.999991993 0.999999207

1.E − 05 1.E − 06 1.E − 07 1.E − 08 1.E − 09 1.E − 10

0.999999638 0.999999256 0.999977639 1.000120637 0.999324411 0.986622386

506

I. Charpentier / Journal of Volcanology and Geothermal Research 175 (2008) 501–508

suffices to modify the printed outputs to u0d, theta0d, r0d, n0d at the end of the calM_D routine. Sensitivity computations performed in Charpentier and Espíndola (2005a,b), may be produced to analyze the Woods model response to change in its boundary conditions within El Chichòn volcanic context. One notices that sensitivity analyses with respect to parameters require the corresponding differentiation of the Plinian code. 4.2. A perfect identical twin experiment Identical twin experiments are data assimilation experiments involving synthetic data created within a model M from known inputs. In the following, the Plinian model is run from inflow condition u# 0 = (300 m/s, 1000 K, 100 m, 0.03) for El Chichón geographical data (1060 m a.s.l.) and a standard atmosphere (tropopause at 16,500 m a.s. l.). Such a simulation generates a complete set of perfect data. Previous adjoint experiments on eruptive columns (Charpentier, 2007) show that observation sets constituted of velocity, temperature and radius data ensure the identification of the inflow conditions. This happens

even an important white Gaussian noise is present on the data. Nevertheless, the larger the noise, the larger the data set to satisfy the white Gaussian noise hypothesis. As presented in the former paper, the identification may occur when only a couple of variable (and sometimes a unique variable) is under observation. The PliniAD objective function computes discrepancies on velocity, temperature and radius data each 500 m between the z0 + 500 m and the tropopause. If the user makes changes to routine calJ, a new differentiation of the code is needed to guarantee the correctness of the adjoint code. The PliniAD package includes a perfect twin experiment that allows to prove the correctness of the optimization process. The latter is driven by the bound-constrained quasi-Newton algorithm (L-BFGS-B library) proposed by Lu et al. (1995), the gradient being computed with the adjoint code. In order to avoid collapsing columns during the iterative process, the space of control (see the driver routine) uses typical ranges for Plinian columns: u0 ∈ (200 m/s, 600 m/s), θ0 ∈ (800 K, 1300 K), r0 ∈ (50 m, 200 m) and u0 ∈ (0.02, 0.07). The twin_experiment call graph is organized as follows:

Fig. 1. Optimal and identified velocities (m/s), radii (m), gas fractions and temperature (K) at the vent, optimal and identified entrainment parameters, and discrepancy results for the 10 identical twin experiments.

I. Charpentier / Journal of Volcanology and Geothermal Research 175 (2008) 501–508

In the package, the twin_experiment routine first generates synthetic observations running the Plinian model routine calM from inputs u# 0 . Then the user chooses an initial iterate, for instance u0 = (250 m/s, 900 K, 130 m, 0.04), and performs the optimization calling the driver routine. The latter was adapted from a generic routine provided by the L-BFGS-B library. One notices that driver is an iterative scheme that successively calls setulb (L-BFGS-B library) to update the iterates, and the adjoint code calM_B (of calM) to calculate the objective function and its gradient. As data are perfect ones, the identification is perfect: at the end of the computation, the discrepancy j is less than 10− 10. This twin experiment requires less than 4 s on an Intel Xeon 3.20 GHz. The computation can be done with the tangent linear code since the dimension of the control variable is small. It requires about twice the CPU time because the tangent linear code has to be called 4 times. Another solution would be to generate the tangent linear code using the “Tangent Mutidirectional Mode” proposed by Tapenade, the trajectory being then computed once only. 4.3. Twin experiments involving noisy data Realistic data are obtained by adding white Gaussian noises (of law N (0, σ)) and by sampling the data by means of observation operators. Such a numerical approach allows to avoid modelling errors and to take under control measurement uncertainties. In this paper, the data assimilation process is applied to the identification of both the inflow conditions and the entrainment parameter. This requires a differentiation with respect to these 5 controls. As column data are not available for the El Chichón eruptions, VDA abilities are studied by means of twin experiments, data being generated running the model M from known inputs. Volcanic and atmospheric parameters are those of the PliniAD data.inc file. In order to mimic the temporal evolution of the column, the model was run from 10 inflow conditions arbitrarily chosen such that the boundary mass flux F0 = β0u0r20 is constant. Boundary temperature (θ0 = 1000 K) and gas fraction (n0 = 0.03) remain constant whereas a decrease in boundary velocity is related to an increase in boundary radius. Thus 10 snapshots of the column were produced using increasing values r0 = 75i (i = 0,…, 9) (and corresponding velocity values) to mimic the erosion of the conduit. Generated velocity, temperature and radius data were then perturbed using increasing white Gaussian noises with standard deviations respectively equal to 3 im/s, iK, and 6 im (i = 0,…, 9). With regard to Charpentier (2007), observation is concerned with the tropospheric part of the column by choosing 40 data at heights zd = z0 + 75 im (i = 1, 40). A peculiar attention was put to verify that the white Gaussian noise hypothesis is satisfied for each of the 3 samples of 40 data. As presented in Fig. 1, the VDA process enables the identification of inflow conditions suggesting that it could be possible to follow the temporal evolution of boundary conditions during the eruption. The last picture presents discrepancy results obtained at the end of the optimization process. Even large noises imply large discrepancies, the identification occurs since, as expected from the Statistics theory, data were pertubated with white Gaussian noises and large enough samples were used. 5. Conclusions Both the modelling and the observation of Plinian columns are subject to errors. On the one hand, 1-dimensional one-phase models may appear to be unsuited for the assimilation of data collected on a 3dimensional turbulent column affected by the wind. On the other

507

hand, data require peculiar treatments to be interpreted in the framework of the modelling. Obviously, complementary theoretical and experimental works will have to be performed before the assimilation of actual data set can be applied. Although no column data are available for the 1982 eruptions, VDA experiments presented in this paper show that such data could have been handled for the identification of the inflow conditions. Performing VDA experiments with this Plinian model coupled to a tephra-fall model (Carey and Sparks, 1986) could be an interesting solution to exploit isopach and isopleth maps of the past eruptions of El Chichón (Espindola et al., 2000). Acknowledgments The author wish to thank the reviewers and J. M. Espíndola (Instituto de Geofísica, UNAM) for fruitful discussions. The help from G. Zenteno (Instituto de Geología, UNAM) is gratefully acknowledged. References Bunge, H.P., Hagelberg, C.R., Travis, B.J., 2003. Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography. Geophys. J. Int. 152, 280–301. Calder, E.S., Sparks, R.S.J., Woods, A.W., 1997. Dynamics of co-ignimbrite plumes generated from pyroclastic flows of Mount St. Helens (7 August 1980). Bull. Volcanol. 58, 432–440. Canul, R.F., Rocha, V.L., Informe geológico de la zona geotermica de “El Chichónal”, Chiapas, Com. Fed. de Electr., Morelia, México, 1981. Carey, S., Sparks, R.S.J., 1986. Quantitative models of the fallout and dispersal of tephra from volcanic eruption columns. Bull. Volcanol. 48, 109–125. Carey, S., Sigurdsson, H., 1986. The 1982 eruptions of El Chichón volcano, Mexico (2): observations and numerical modelling of tephra-fall distribution. Bull. Volcanol. 48, 127–142. Charpentier, I., 2000. The MesODiF package for gradients computations with the atmospheric model Meso-NH. Environ. Model. Softw. 15, 533–538. Charpentier, I., 2007. Adjoint modelling experiments on eruptive columns. Geophys. J. Int. 169, 1356–1365. Charpentier, I., Espíndola, J.M., 2005a. Local sensitivity analysis of a numerical model of volcanic Plinian columns through automatic differentiation. Math. Geol. 37, 95–113. Charpentier, I., Espíndola, J.M., 2005b. A study of the entrainment function in models of Plinian columns: characteristics and calibration. Geophys. J. Int. 160, 1123–1130. Chavent, G., 1975. History matching by use of optimal control theory. Soc. Pet. Eng. J. 15, 74–86. Courtier, P., Talagrand, O., 1987. Variational assimilation of meteorological observations with the adjoint equation — part II. Numerical Results. Quart. J. Roy. Meteor. Soc., vol. 118, pp. 649–672. Dobran, F., Neri, A., Macedonio, G., 1993. Numerical simulation of collapsing volcanic columns. J. Geophys. Res. 98, 4231–4259. Espindola, J.M., Macias, J.L., Tilling, R.I., Sheridan, M.F., 2000. Volcanic history of El Chichón Volcano (Chiapas, Mexico) during the Holocene, and its impact on human activity. Bull. Volcanol. 62, 90–104. Griewank, A., 2000. Evaluating derivatives: principles and techniques of algorithmic differentiation. Frontiers in Appl. Math, vol. 19. SIAM, Pennsilvany. Hascoët, L., Pascual, V., 2004. TAPENADE 2.1 user's guide. INRIA technical report RT0300. Herzog, M., Graf, H.F., Textor, C., Oberhuber, J.M., 1998. The effect of phase changes of water on the development of volcanic plumes. J. Volcanol. Geotherm. Res. 87, 55–74. Hoblitt, R.P., 1986. Observations of the Eruptions of July 22 and August 7, 1980, at Mount St. Helens. Professional Paper, vol. 1335. U.S. Geological Survey, Washington, pp. 1–44. Ide, K., Courtier, P., Ghil, M., Lorenc, A.C., 1997. Unified notation for data assimilation: operational, sequential and variational. J. Met. Soc. Japan 75, 181–189. Ishimine, Y., 2006. Sensitivity of the dynamics of volcanic eruption columns to their shape. Bull. Volcanol. 68, 516–537. Kieffer, S., 1981. Blast dynamics at Mt. St. Helens on 18 May 1980. Nature 291, 568–570. Kieffer, S., Sturtevant, B., 1981. Laboratory studies of volcanic jets, 1984. J. Geophys. Res. 89, 8253–8268. Le Dimet, F.-X., Talagrand, O., 1986. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. Tellus 38, 97–110. Lewis, J., Derber, J., 1985. The use of the adjoint equations to solve a variational adjustment problem with advective constraints. Tellus 37, 309–327. Lions, J.-L., 1971. Optimal control of systems governed by partial differential equations. Springer-Verlag, Berlin. Lu, P., Nocedal, J., Zhu, C., 1995. A limited memory algorithm for bound constrained optimization. SIAM J. Sci. Comput. 16, 1190–1208. Matson, M., 1984. The 1982 El Chichón Volcano eruptions — a satellite perspective. J. Volcanol. Geotherm. Res. 23, 1–10. Morgensten, J., 1985. How to compute fast a function and all its derivatives, a variation on the theorem of Baur–Strassen. SIGACT News 16, 60–62. Morton, B.R., Taylor, G., Turner, J.S., 1956. Turbulent gravitational convection from maintained and instantaneous sources. Proc. R. Soc. Lond., A 234, 1–23.

508

I. Charpentier / Journal of Volcanology and Geothermal Research 175 (2008) 501–508

Nairn, I.A., 1976. Atmospherical shock waves and condensation clouds from Ngauruhoe explosive eruptions. Nature 259, 190–192. National Earth Satellite Service (Ness), 1982. El Chichón volcanic eruptions (WAB467, 10 minutes silent film). Walther Boham Co., Washington, DC. Oberhuber, J.M., Herzog, M., Graf, H.F., Schwanke, K., 1998. Volcanic plume simulation on large scale. J. Volcanol. Geotherm. Res. 87, 29–53. Pironneau, O., 1983. Optimal shape design of elliptic systems. Lecture Notes in cpm. phys. Springer. Prandtl, L., 1949. The essentials of fluid dynamics. Blackie & Son Ltd, Glasgow. Robock, A., 2000. Volcanic eruptions and climate. Rev. Geophys. 38, 191–219. Roche, J.R., 1997. Gradient of the discretized energy method and discretized continuous gradient in electromagnetic shaping. Appl. Math. Comput. Sci. 7, 545–565. Settle, M., 1978. Volcanic eruption clouds and the thermal power output of explosive eruptions. J. Geophys. Res. 3, 309–324. Sigurdsson, H., Carey, S.N., Espindola, J.M., 1984. The 1982 eruptions of El Chichón volcano, Mexico: stratigraphy of pyroclastic deposits. J. Volcanol. Geotherm. Res. 23, 11–38. Sirkes, Z., Tziperman, E., 1997. Finite difference of adjoint or adjoint of finite difference? Mon. Weather Rev. 125, 3373–3378. Sparks, R.S.J., 1986. The dimension and dynamics of volcanic eruption columns. J. Volcanol. Geotherm. Res. 48, 13–15.

Sparks, R.S.J., Bursik, M.I., Carey, S.N., Gilbert, J.S., Glaze, L.S., Sigurdsson, H., Woods, A.W., 1997. Volcanic plumes. Wiley J. & Sons, Chichester. Talagrand, O., Courtier, P., 1987. Variational assimilation of meteorological observations with the adjoint vorticity equation. I - Theory. Quart. J. Roy. Meteor. Soc. 113, 1311–1328. Textor, C., Graf, H., Longo, A., Neri, A., Ongaro, T.E., Papale, P., Timmreck, C., Ernst, G.G.J., 2005. Numerical simulation of explosive volcanic eruptions from the conduit flow to global atmospheric scales. Ann. Geophys. 817–842. Valentine, G.A., Wolhetz, K.H., 1989. Numerical models of Plinian eruption columns and pyroclastic flows. J. Geophys. Res. 94, 1867–1887. Wilson, L., 1976. Explosive volcanic eruptions III. Plinian eruption columns. Geophys. J. R. Astron. Soc. 45, 543–556. Wilson, L., Sparks, R.S.J., Huang, T.C., Watkins, N.D., 1978. The control of volcanic column heights by eruption energetics and dynamics. J. Geophys. Res. 83 (B4), 1829–1836. Woods, A.W., 1988. The fluid dynamics and thermodynamics of eruption columns. Bull. Volcanol. 50, 169–193.