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.