Finite element models of volcano deformational systems having structural complexity

Resumen Abstract Índice Conclusiones

Ronchin, Erika

2016-A
Descargar PDF

Resumen

«FINITE ELEMENT MODELS OF VOLCANO DEFORMATIONAL SYSTEMS HAVING STRUCTURAL COMPLEXITY»

Erika Ronchin

Department of Geodynamics and Geophysics

Universitad de Barcelona

Barcelona 2015

Resumen

El objetivo principal de este trabajo es la construcción de modelos de elementos finitos (FEMs) 3-D con complejidades estructurales con el fin de simular sistemas volcánicos de manera más realista. Como ejemplo de aplicación se ha escogido la caldera de Rabaul, un sistema volcánico cuya dinámica no se comprende por completo. Invirtiendo los datos de InSAR recogidos durante los años 2007-2010, investigamos las fuentes de desplazamiento de la superficie y proporcionamos claves de relevancia sobre el sistema magmático superficial real.

Incluyendo características realistas, como la topografía y heterogeneidades mecánicas, usamos las informaciones geofísicas y geológicas para construir modelos de FEMs complejos en 3D. En última instancia, proporcionamos una estrategia para llevar a cabo una inversión lineal basada en una matriz de fuentes que nos permite identificar una distribución de flujo de fluido a través de un volumen de posibles fuentes responsables de los cambios de presión en el medio según lo dictado por los datos, sin imponer a priori una forma de fuente específica y su profundidad.

Los resultados permiten generar imágenes de la forma compleja de la fuente que da lugar a la deformación, en el espacio y en el tiempo, sin tener que utilizar ninguna fuente con una forma excesivamente simplificada a priori. Esto lleva el modelado de fuentes un paso adelante hacia modelos más realistas.

En el caso de Rabaul, la aplicación de la metodología discutida anteriormente, muestra un sistema magmático superficial formado por dos lóbulos interconectados localizados bajo la caldera y en posiciones diametralmente opuestas. La interconexión y la distribución espacial de las fuentes encuentran correspondencia en la petrología de los productos descritos en literatura y en la dinámica de las erupciones que caracterizan la caldera. Los resultados obtenidos mediante la aplicación del método son satisfactorios y demuestran que la inversión lineal basada en la matriz de fuentes de FE propuesta puede ser considerada adecuada para generar modelos de sistemas magmáticos. Se puede aplicar fácilmente a cualquier volcán, ya que tiene en cuenta la deformación del edificio sin tener que especificar la forma de la fuente de deformación antes de la inversión.

Abstract

Abstract of PhD thesis: «FINITE ELEMENT MODELS OF VOLCANO DEFORMATIONAL SYSTEMS HAVING STRUCTURAL COMPLEXITY»

Erika Ronchin

Department of Geodynamics and Geophysics

Universitad de Barcelona

Barcelona 2015

Abstract:

The main focus of this work is to build 3-D FEM models with structural complexities in order to simulate volcanic systems in a more realistic way. We use Rabaul as an example to show the application of the methods and strategies proposed to an active volcano.

Rabaul caldera is a volcanic system whose dynamics still need to be understood to effectively predict the behavior of future eruptions. In comparison to the simplified analytic models used so far, more realistic models, such as Finite Elements Models (FEMs), are needed to more accurately explain recent deformation and understand the magmatic system. By inverting InSAR data collected between 2007 and 2010 (using linear inversions based on FEMs), we investigate the sources of surface displacement and provide insights about the actual shallow magmatic system.

FEMs are numerical models that let us include realistic features such as topography and mechanical heterogeneities. We provide strategies to use geophysical and geological information to build complex 3-D parts and assemble them into 3-D models. We then compare the effects of different material properties configurations and of different source shapes on the deformational signal and on the strength source estimates (fluid flux or pressure).

Ultimately, we provide a strategy for performing a linear inversion based on an array of FEM sources that allows us to identify a distribution of flux of fluid (or change in pressure) over a volume, without imposing an a-priori source shape and depth. We use Rabaul as an example to show the 3-D model’s validity and applicability to active volcanic areas. The methodology is based on generating a library of forward numerical displacement solutions, where each entry is the displacement generated by injecting a mass of fluid of known density and bulk modulus into a source of the array. The sources are simulated as fluid-filled cavities that can accept a specified flux of magma. As the array of sources is an intrinsic geometric aspect of all forward models and the sources are activated one at a time, the domain only needs to be discretized once. This strategy precludes the need for remeshing for each activated source and greatly reduces computational requirements. By using an array of sources, we are not investigating the geometric and pressure parameters of a simplified, unique source with a regular shape. Instead, we are investigating a distribution of flux of fluids over a volume of potential sources responsible for the pressure changes in the medium as dictated by the data. The results allow us to image the complex shape of the deformation source without having to use any a-priori or simplified sources. This takes source modeling a step towards more realistic source models.

The application of the methodology to Rabaul shows a shallow magmatic system under the caldera made of two interconnected lobes located at the two opposite sides of the caldera. These lobes are suggested to be the feeding reservoirs of the ongoing Tavuvur volcano eruption, on the eastern side, and of the past Vulcan volcano eruptions, on the western side. The interconnection and spatial distribution of sources find correspondence in the petrography of the products described in literature and in the dynamics of the single and twin eruptions that characterize the caldera.

The good results obtained from the application of the method show that the proposed linear inversion based on the FEM array of sources can be considered suitable for generating models of the magmatic system. It can be easily applied to any volcano, because it accounts for volcano deformation without having to specify the shape of the deformation source prior to inversion.

Índice

Table of contents

Abstract and Resumen . . . . . . Pg.1

Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

Resumen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1 Introduction . . . . . . Pg.4

1.1 Motivation and objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.1.1 Rabaul caldera . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.1.2 Volcano models: state of the art . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.1.3 Rabaul models: state of the art . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.1.4 Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.2 Approach to the problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1.3 Structure of the dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2 Rabaul geology, volcanic activity, and data . . . . . . Pg.13

2.1 Tectonic setting of New Britain and Gazelle Peninsula . . . . . . . . . . . . . . . . . 13

2.1.1 Regional tectonic evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.1.2 Stratigraphy and main structures of New Britain Island and Gazelle peninsula 17

2.2 Rabaul caldera geological setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.2.1 Volcanoes pre- and post- caldera formation . . . . . . . . . . . . . . . . . . . 26

2.2.2 Evolution of the magmatic system and actual magma chamber . . . . . . . . 28

2.2.3 Inside the caldera: bulge, circular seismicity, and caldera structures . . . . . . 30

2.2.4 A wide picture of the study area through tomography and material velocities inferred from the tomography . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.3 Rabaul caldera eruptive activity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

2.3.1 Historical eruptions up to the 1980s’ crisis and 1994 twin eruption . . . . . . 34

2.3.2 Recent volcanic activity: period from February 27, 2007 to December 8, 2010 . . . 35

2.4 Rabaul caldera data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

2.4.1 Relief data (topography and bathymetry) . . . . . . . . . . . . . . . . . . . . 38

2.4.2 ALOS PALSAR and PS-InSAR basic concepts . . . . . . . . . . . . . . . . . 39

ALOS PALSAR basic concepts . . . . . . . . . . . . . . . . . . . . . . . . . . 40

PS-InSAR technique basic theoretical concepts . . . . . . . . . . . . . . . . . 44

2.4.3 Rabaul ALOS-PALSAR data: PS images, mean velocity, and velocity standard deviation . . . . 45

Insights from the temporal series . . . . . . . . . . . . . . . . . . . . . . . . . 46

Insights from cumulative displacements . . . . . . . . . . . . . . . . . . . . . 49

Insights from mean velocity and its STD . . . . . . . . . . . . . . . . . . . . . 50

3 Methods and procedures . . . . . . Pg.53

3.1 Analytical models and finite element models (FEMs) . . . . . . . . . . . . . . . . . . 53

3.1.1 Analytical models for volcanoes: point source, finite spherical pressure source, and others . . . . . . 53

3.1.2 Finite element models (FEMs) . . . . . . . . . . . . . . . . . . . . . . . . . . 56

3.2 Modeling protocol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3.2.1 Define purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3.2.2 Conceptual model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

3.2.3 FEM model design and configuration . . . . . . . . . . . . . . . . . . . . . . 61

3.2.4 Calibration (inverse analysis) . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

3.2.5 Verification and Post-audit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

3.3 Building strategy for 3-D parts in Abaqus . . . . . . . . . . . . . . . . . . . . . . . . 63

3.3.1 3-D modeling strategy using Python script to implement Abaqus CAE . . . 63

3.4 Effective material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

3.4.1 Poisson’s ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

3.4.2 Young’s modulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

3.4.3 Dynamic elastic moduli ndy and Edy, and density . . . . . . . . . . . . . . . . 70

3.4.4 Upscaling of dynamic moduli to more appropriate properties . . . . . . . . . 70

3.4.5 Effective magma bulk modulus, b . . . . . . . . . . . . . . . . . . . . . . . . . 73

3.5 Mesh model validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

3.5.1 The importance of mesh validation . . . . . . . . . . . . . . . . . . . . . . . . 75

3.5.2 Estimating discrepancies between the FEM and reference model, and considering an appropriate level of mesh resolution errors . . . . . . . . . . . . . . . 76

3.5.3 About the importance of choosing the appropriate reference model and nodes for mesh testing . . . . 77

3.6 Weighted Quadtree-based algorithm for data reduction . . . . . . . . . . . . . . . . 81

3.6.1 Introduction to sub-sampling procedures . . . . . . . . . . . . . . . . . . . . . 81

3.6.2 Modified variance-quadtree algorithm (VQT) . . . . . . . . . . . . . . . . . . 83

3.7 Linear inverse theory and strategies . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

3.7.1 Linear least squares solution of the inverse problem . . . . . . . . . . . . . . . 87

3.7.2 Improvements of model parameters estimation . . . . . . . . . . . . . . . . . . 89

3.7.3 The finite difference laplacian operator (3D discrete smoothing operator of regular grid sources with zero Dirichlet boundary condition) . . . . . . . . 92

3.8 FE array of sources at the base of the Green’s function matrix generation . . . . . . 95

3.8.1 Cubic source generation: Fluid elements and aspects of the cavity construction algorithm . . . . . . 95

3.8.2 Array of source-elements and verification of the cubic pressure sources with McTigue . . . . . . 99

4 Results of methods applied to Rabaul caldera . . . . . . Pg.104

4.1 Rabaul geologic parts identification and 3-D construction in Abaqus CAE environment104

4.1.1 Model mantle (50-30 km depth) and lower crust (30-8 km depth) . . . . . . . 105

4.1.2 Model parts of the upper crust (8-0 km depth) . . . . . . . . . . . . . . . . . 106

Intra-caldera fill . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

Dikes swarm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

Baining Mountains Complex . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

Extra-caldera sediments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

Rabaul magma chamber and magma bulk modulus, b . . . . . . . . . . . . . 114

4.2 Results of the Rabaul 3-D forward models . . . . . . . . . . . . . . . . . . . . . . . 119

4.2.1 Assembling the parts in a 3-D model with topography . . . . . . . . . . . . . 119

4.2.2 Boundary conditions and mesh validation . . . . . . . . . . . . . . . . . . . . 120

4.2.3 Results of 3-D Rabaul Abaqus models . . . . . . . . . . . . . . . . . . . . . . 124

Preliminary study of the topographic effects on surface deformation related to changes in material properties. . . . . . 126

Effects of the Poisson’s ratio and Young’s modulus . . . . . . . . . . . . . . 127

Influence of the geologic bodies-single parts . . . . . . . . . . . . . . . . . . . 132

Combined parts effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

Magma chamber shape effects . . . . . . . . . . . . . . . . . . . . . . . . . . 142

Combined effects of magma chamber shape and soft caldera infill . . . . . . 145

All results summarized . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147

4.3 Rabaul PS-InSAR quadtree reduction . . . . . . . . . . . . . . . . . . . . . . . . . . 151

4.4 Results of Rabaul inverse models using a-priori sources . . . . . . . . . . . . . . . . . 154

4.5 Results of Rabaul inverse model using a 3-D source array . . . . . . . . . . . . . . . 160

5 Discussion . . . . . . Pg.171

5.1 Geometry of the models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171

5.2 On the importance of choosing the material properties . . . . . . . . . . . . . . . . . 172

5.3 Effects of Rabaul material properties on the displacements signals . . . . . . . . . . 173

5.4 Model assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175

5.5 Influence of material properties and magma chamber geometries on deformation predictions and pressure estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175

5.6 Search of magma withdrawal through the inversion based on FE array of sources . . . . . . . . 180

6 Conclusions, recommendations and future works . . . . . . Pg.187

6.1 Conclusions and recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187

6.2 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188

Bibliography . . . . . . Pg. 190

Appendix A CAESPLINE.py . . . . . . Pg.211

Appendix B IDL procedure to compile the Laplacian operator . . . . . . Pg.214

Conclusiones

Conclusions, recommendations and future works

7.1 Conclusions and recommendations

In this work, methods and strategies are developed to build a 3-D

FEM with structural complexities and integrate it into a linear

least squares damped analysis. The methods are applicable to any

volcanic area. Here, all methods and models are applied to Rabaul

caldera with the main goal of imaging the shallow magma chamber,

intended here as a volume of the model through which fluids can

migrate. The application to a real case motivated us to study

many of the aspects of the model in detail, coming up with

general considerations about modeling natural volcanic areas in

more realistic ways.

The topography needs to be taken into account, even in the case

of non-prominent features, whenever steep slopes are present in

the vicinity of the sources, such as in the case of Rabaul

caldera. In case very weak material is constituting the upper

part of the model (e.g. the saturated material of the Rabaul

caldera infill) and no in-situ elastic properties estimation is

available, a sensitivity test to the Poisson’s ratio is

recommended. The fact that a sharp transition of material

properties closer to the source impacts the deformation more

highlights the importance of modeling the medium with right

material properties close to the source. This means that when the

source is deeper, we also need to make sure to model the deeper

material well. Furthermore, a magma chamber inferred from the

tomography and seismic data does not predict the displacements

better than a spherical one with the same volume. This is due to

the fact that, although using a complex geometry, we have the

limitation that we can only apply one single pressure uniformly

over the walls of the cavity. Due to this limitation and to the

ambiguity of the sources of low tomographic velocities, it is

fairly difficult to achieve a good model of the magma chamber

modelled as a cavity based on tomographic images. Nevertheless,

the tomography provides a great guide to define a volume of

investigation for the array of sources at the base of the method

proposed in this dissertation.

We propose a method to implement a linear least squares inversion

with FEM for the estimation of an arbitrary distribution of mass

change among a volume of the FEM. This leads to the possibility

of imaging the magma chamber through the inversion of geodetic

data, thus providing more insights about the magmatic system with

respect to the over-simplified a-priori cavity used to simulate

the magma chamber. The method pushes FEM-based inverse models a

step forward with the result of being able to image the magmatic

system in a more realistic way through simple linear inversion

strategies.

By applying the strategies and methods to the active Rabaul

caldera, we are able to image the extension of the magma chamber

under the caldera. This could have implications for future

evaluation of the possibility of catastrophic collapses. In fact,

the imaged volume, included in the area defined by the shallow

seismicity, looks to be sufficient to explain the observed

deformation, with no need of a bigger magma chamber as it would

be suggested by the laterally extended low velocity contrast of

the tomographic images. This has implications for hazard

assessment in the sense that the dimension of a possible

worst-case scenario could be estimated.

Our study confirms the presence of two main sources at opposite

sides of the caldera hypothesized by previous analytical models,

and images the interconnection between these two sources. Thus,

the Rabaul shallow magmatic system seems to act like a unique

reservoir instead of different disconnected magma reservoirs.

This implies a high mobility of magma between different parts of

the system. Thus, as shown by the deformational period studied, a

big eruption at one of the centers, Tavuvur in this case, can

trigger a long term deflation of the entire system, Vulcan batch

included. The imaged connection could have implications for the

dynamics of future twin eruptions. Our study may thus help to

understand the dynamics of this kind of eruptions and forecast

them.

7.2 Future works

Although this study images the magma chamber under Rabaul

caldera, we have to remember that border effects have been

observed on the estimates of the inversion, especially on the

northern side of the FE source array. In future works, a bigger

array of sources could be used in order to avoid border effects,

to be able to investigate the full low velocity zone shown by the

tomography of Finalyson et al. (2003), and to extend the study to

deeper depths in order to investigate the effects of deeper

reservoirs on the surface displacement signal. Furthermore,

although not much of the hydrothermal system is known until now (Johnson et al., 2010)

, future studies will perhaps allow us to extend the model to

shallower levels in order to investigate the effects on the

hydrothermal system of the surface deformations and better

predict the area north of Tavurvur.

The estimated distribution of pressures looks reasonable,

encouraging us to use the model in future studies for the

understanding of stress changes and failure mechanisms of the

area. A more realistic model in terms of stress distribution

should also include geostatic conditions.

In this study, the domain in which the sources are included is

modeled using elastic materials inferred from seismic velocities

and geologic observations. Future studies could investigate

alternative rheologies (e.g. visco-elasticity) to better describe

the domain rheology.

At this point, the study is limited to the inversion only of 1-D

ascending InSAR data. For future studies, we hope that we will

have access to a 3-D description of the signal, and thus hope to

integrate the InSAR data with GPS data in order to have a more

robust inversion.

As our models predict mass movements through the magmatic system,

in future studies we could also calculate the related gravity

anomalies and introduce the gravity anomalies measured at the

surface into the inverse scheme. By considering these anomalies,

we could detect processes other than mass flux that contribute to

the surface deformation without requiring a change of mass (e.g.

thermoelastic expansion).

Future studies could also combine the inversion of surface

deformation and gravity anomalies with the monitored seismicity

to track magma migration though the system.

By referring to the mapped flux of magma when interpreting the

surface deformation, we suggest that the deformation is only due

to the movement of magma within the plumbing system. Besides the

movement of magma, a wide range of magmatic processes could be

responsible for the pressure variations (or volume variations)

that generate surface deformation (e.g. crystallization,

degassing, expansion of the hydrothermal system, remelting,

etc.). As the change of pressure inside the elements and the

change of volume of the elements are solutions of the forward

unity flux models and they are linearly correlated to the flux

solutions, we can provide maps of pressure change (Fig.5.6.1)

and volume change. By comparing these maps with values of

pressure change (or volume change) due to other magmatic

processes, it will be possible in the future to balance the

effects of different processes and provide a more complete and

realistic overview of the processes responsible for the observed

deformation.

At small time intervals the response of the system can be

considered elastic. An analysis of short temporal intervals could

therefore take full advantage of the model here proposed. With a

temporal analysis of surface displacements based on one single

FEM, we can answer questions about how the magma is moving in the

system, if it is being accommodated by the entire system, or if

it is accumulating on one side, building up stresses on a

specific area. Answers to these questions are fundamental to

understand the running dynamics in the system and to forecast the

next activity.

Lastly, the 3-D model proposed for Rabaul caldera could be used

in new tomographic studies to better find the anomalies of

velocity.