Imaging the complex geometry of a magma reservoir using FEM-based linear inverse modeling of InSAR data: application to Rabaul Caldera, Papua New Guinea

Erika Ronchin,1,2 Timothy Masterlark,3 John Dawson,4 Steve Saunders5 and Joan Martı̀ Molist1 1Group of Volcanology of Barcelona, Institute of Earth Sciences Jaume Almera, ICTJA-CSIC, Lluis Solé i Sabaris s/n, Barcelona E-08028, Spain. E-mail: erikaronchin@gmail.com 2Mineralogy Petrology Tectonics Department of Earth Sciences, Uppsala University, Villavägen 16, SE-75236 Uppsala, Sweden 3Department of Geology and Geological Engineering, South Dakota School of Mines, Rapid City, SD 57701, USA 4Geodesy and Seismic Monitoring Group, Geoscience Australia, Cnr Jerrabomberra Ave and Hindmarsh Drive, Symonston, ACT 2609, Australia 5Rabaul Volcano Observatory, P.O. Box 386, Rabaul, East New Britain Province, Papua New Guinea

Imaging a magma reservoir using FEMs and InSAR 1747

I N T RO D U C T I O N
A significant challenge in assessments of volcano hazards is the accurate imaging of deformation sources, which represent magma migration and storage in both space and time. Our ability to predict future eruptions also relies on the knowledge of which parts of a reservoir are active versus those that could reactivate in the future due to magma migration through connected parts of a composite reservoir. Using customary simplistic a-priori deformation source geometries to represent these aspects strongly limits our explorations. The goals of this study, therefore, are (1) the introduction and testing of an innovative inversion that is based on finite-element method (FEM) models and that does not require an a-priori source geometry and (2) the application of this approach to the imaging of the extended shallow magmatic system responsible for the broad and long-term subsidence of Rabaul caldera. Mathematical models provide a link between the sources at depth (e.g. the deforming magmatic reservoir) and the observed system response at the surface. Parameters that describe the source (e.g. strength and geometric parameters) can be estimated, through inverse modeling, by comparing the observed natural system response (e.g. surface displacement) to the one predicted by a mathematical forward model of the source. The reliability of these models is customarily based on the goodness of predictions for a given set of available information, such as an interferometric synthetic aperture radar (InSAR) image.
Models with more complex source geometries generally allow the estimation of more parameters. This makes possible a more detailed investigation of real sources; however, at higher computational cost than simplified models. Therefore, subsurface processes of very complex volcanic systems are still interpreted on the basis of oversimplified sources that may ignore the complete assembly of available information. A consequence of this approximation is the estimation of a uniform, maximum pressurization of the magma reservoir which is incompatible with the rock strength. The pressurization can be reduced to more realistic values by claiming the effects of the inelasticity of the medium (e.g. Fernández et al. 2001;Trasatti et al. 2003;Hickey et al. 2015).
FEM models, referred to in the following as FEMs, allow us to take into account many realistic aspects of volcanic systems (e.g. topography and mechanical heterogeneities) which are difficult to represent with analytical models. However, the potential for the application of FEM models is restricted because of the use of computationally convenient source simplifications. Common simplifications of source geometry imply the use of a-priori and/or simplified source shapes. Based on geophysical studies, such as tomography, and geological field evidence, it is well known that magmatic reservoirs are generally geometrically complex (e.g. Smith et al. 2009;Burchardt et al. 2010Burchardt et al. , 2013Jaxybulatov et al. 2014;Huang et al. 2015). This means that unrealistic aspects of the volcanic systems are introduced in the analysis by using models with simplified sources.
Using the principle of superposition of point sources embedded in elastic homogeneous half-spaces, a few studies attempted to simulate complex sources having arbitrary geometries (Vasco et al. 2002;Masterlark & Lu 2004). However, their results were limited by the simplification introduced when representing the complex heterogeneous nature of the volcanic edifice through an elastic homogeneous half-space. The goal of these studies was the characterization of the unknown source without using a-priori geometry, solving for the source strengths of an amorphous cluster of deformation sources distributed over a grid (e.g. Mossop & Segall 1999;Vasco et al. 2002;Masterlark & Lu 2004;Camacho et al. 2011). This approach involves an array of simple sources, in which the source strengths are linearly related to the deformation via elastic equations. This allows the application of linear inverse methods.
The assumption of linearity is common in volcano models, partly because the data may not constrain more complex rheology (Trasatti et al. 2008). Also, modeling inelastic behaviour requires much more computation.
Such a linear inversion based on an array of sources gives a geodetic image (Vasco et al. 2002) made of a distribution of source strengths that allows the identification of the position and the arbitrary shape of the reservoir. Presumably, this leads to a more realistic interpretation of deformational sources having complex geometry. This approach is generally applied to analytical models, for which the solutions used to build the Green's functions matrix are quickly and easily computed. The strategy is computationally challenging if applied to an inverse scheme based on FEMs, due to the fact that a new mesh of the domain needs to be generated for each source. Trasatti et al. (2008) overcame this problem by loading the faces of a cubic element of a meshed domain with three dipoles and three double couple forces, assembled in a stress tensor. However, their aim was to find depth, strength and ellipsoidal shape of a single pressurized source through the estimation of six elements of the stress tensor. If applied to the estimation of a grid of stress tensors, this source formulation would still generate a very expensive computational problem due to the high number of parameters to be estimated. Also, it needs to be justified by the quantity and the quality of available data sets (Trasatti et al. 2008). As a result, their investigation was limited to a single-source tensor.
We provide a computationally inexpensive strategy for the linear regularized least-squares inversion of InSAR data based on the Green's functions of an extended grid of FEM sources. The goal of the proposed inversion is to study the withdrawal of magma from a reservoir in order to define the amorphous shape and location of the unknown source in a volume previously constrained by the seismicity, geology and tomography. We apply our methodology to InSAR data of Rabaul caldera, Papua New Guinea, which recorded a broad long-term subsidence centred at the caldera between 2007 February and 2010 December. The volume within the FEMs that contains the magma reservoir of unknown shape is modeled with a grid of cubic sources, pressurized by the injection of magma mass, whose strength, implemented here as a mass flux of magma and resulting pressure change, is estimated during the inversion. This provides, through a source strength distribution, a more realistic image of the magma reservoir responsible for recent single and twin eruptions at the caldera.

R A B AU L C A L D E R A D E F O R M AT I O N
Rabaul is a historically active collapse caldera with an elliptic shape having major and minor axes of 14 and 9 km, respectively, and outward dipping seismicity due to the elliptical fault zone located at the centre of the caldera (Fig. 1a). The volcanic activity is characterized by single eruptions (Tavurvur) at the eastern side of the caldera or twin eruptions (Tavurvur-Vulcan) at opposite sides of the caldera.
Analytical models of Rabaul caldera displacements, which occurred during the deformational crisis of 1983-1985 and the period prior to the twin eruptions of 1994 (which culminated in a plinian eruption of Vulcan), estimated a shallow source (between 1.2 and 2 km depth) southeast of Matupit Island, under Greet Harbor (McKee et al. 1984, 1989Archbold et al. 1988), and a deeper source near Vulcan at 3 km depth (McKee et al. 1984). FEMs have been used to explain the deformation prior to the 1994 twin eruptions with the interaction between an expanding reservoir and the elliptical faults (De Natale & Pingue 1993;De Natale et al. 1997;Ronchin et al. 2013) or the migration of fluids into the elliptical faults (Saunders 2001). To the authors' knowledge, no FEMs have been constructed to understand the recent deformation and the actual shallow magmatic system that feeds frequent and important eruptions at Rabaul caldera.
In 2006 October, Tavurvur erupted >0.2 km 3 of andesite with a subplinian eruption (Saunders et al. 2007). A period of quiescence followed, with no surface deformation. From 2007 February, the activity of Rabaul caldera area was characterized by almost continuous Tavurvur vulcanian eruptions and by the deformation modeled in this work: a general long-term subsidence of a broad area centred at the caldera and detected by InSAR data between 2007 February and 2010 December. During this period, no changes of products composition were observed (Rabaul Volcano Observatory (RVO) reports in

Rabaul InSAR data
We employed 20 PALSAR images acquired by the ALOS satellite of the Japanese Space Agency (JAXA) between 2007 February 27 and 2010 December 8 (information summarized in Table 1). We processed the ALOS raw data using the ROI PAC software developed by JPL/Caltech (Rosen et al. 2004). We formed the interferograms using the Doris software developed by the Delft University of Technology with respect to a master image (see Supporting Information, Fig. S3), chosen to minimize perpendicular and temporal baselines. Unlike conventional InSAR analysis, no azimuth and range spectral filtering was undertaken to maintain the resolution ). The topographic signal was estimated using the 3 arcsec digital elevation model from the Shuttle Radar Topography Mission (SRTM, Farr et al. 2007). A time-series analysis of the interferograms was then undertaken using the persistent scatterers (PS) InSAR approach implemented in the StaMPS software ) (percent rand and clap win set to 20 and 32, respectively). The PSInSAR technique has been used to overcome many of the problems associated with decorrelation observed in conventional InSAR techniques by identifying pixels which remain coherent (phase stable) over long periods of time . A major advantage of the PSInSAR method for studies of volcanoes is that no a-priori functional model of the temporal deformation is required, since the method relies on an assumption that any deformation signal is spatially correlated. Deformation was estimated simultaneously with nuisance parameters, including atmospheric and residual topography signals, where the signal components were separated by exploiting their temporal and spatial characteristics. Phase unwrapping was performed using the 3-D approach of  (for time-series, see Supporting Information, Fig. S3). Line of sight (LOS) mean velocities were estimated with respect to a manually selected area of pixels outside the area of active deformation. The uncertainty (1σ ) of the velocities was estimated as ±1 cm yr −1 and determined from an analysis of the individual time slices of pixels outside the area of deformation.
Because one goal of this study is the understanding of the extended shallow magmatic system responsible for the broad and long-term subsidence, we investigate the LOS mean velocities over the entire time period of SAR acquisition (Fig. 1a), assuming that this period is long enough to record the response of the entire magma reservoir to a long, almost continuous magma withdrawal due to the Tavurvur eruptions.
A wide area, including the caldera and an area west of the caldera, shows relative negative LOS velocities, with highest values of −60 mm yr −1 between the southwestern tip of Matupit Island and Vulcan's northeastern side (Fig. 1a). Positive LOS displacements are recorded in the eastern areas outside the caldera border. Standard deviations (SD) associated with the mean LOS velocities show the deviation from a linear steady subsidence over the study period ( Fig. 1c). High SD for Matupit Island identify strong nonlinear motion of PS, possibly related to the most shallow processes (e.g. impulses of magma ascent to the surface or hydrothermal processes). The highest SD, along the southwestern side of the old airport strip, are associated with a localized area of anomalous positive LOS velocities suspected to be related to hydrological processes, since this is a drainage area.
The unwrapped InSAR image contains approximately 22 × 10 3 PS (Fig. 1a). In order to accommodate computational requirements of the modeling, we reduced the number of pixels using a variance quadtree algorithm (Samet 1990;Lohman & Simons 2005;Minasny et al. 2007) that spatially subsamples the interferogram (see Supporting Information). With this procedure, we reduced the number of data by about two orders of magnitude (Table 1)   A profile through the InSAR data ( Fig. 1e) shows the most important features of the signal: a broad region of negative LOS velocities outside the caldera, a high-velocity gradient above the Vulcan edifice, loss of data in the area occupied by water and positive LOS velocities above Turanguna and the southern part of Tavurvur.
Neither the seismicity recorded (RVO reports in Smithsonian Institution 2013), nor the distribution of InSAR velocities, suggests movement along the outward dipping elliptical faults during the study period. This encourages us to investigate an extended subsurface magma withdrawal from the shallow reservoir as a possible source of deformation.

M E T H O D : L I N E A R I N V E R S I O N B A S E D O N F E M s , U S I N G R A B AU L C A L D E R A A S A N E X A M P L E
Inspired by the source cluster method consisting of a relatively dense 3-D array of analytical point sources (Vasco et al. 2002;Masterlark & Lu 2004;Battaglia & Vasco 2006;Camacho et al. 2007Camacho et al. , 2011, we model a volume of the Rabaul FEM model domain with a prism of cubic hexahedral elements (Abaqus 2009) that serves to generate a regular 3-D array of elementary sources (Fig. 2b). These are generated one at a time by removing elements of the prism and pressurizing the resultant cavities by injecting a fluid flux. The collection of all source solutions generates a library of pre-computed elementary solutions on the basis of the Green's functions matrix used to solve the inverse problem. The fluid flux and related pressurization in each cavity are estimated through linear inversion. This way, rather than imposing an a-priori source geometry, we can assume that the source of deformation is represented by a subset of the array of sources with an amorphous shape dictated by the data. The use of FEMs for this process allows us to simultaneously account for the arbitrary geometric configurations and the distributions of elastic properties.

FEMs configuration and deformation sources
New FEMs for Rabaul caldera (Fig. 2) are built by implementing a model having heterogeneous distribution of elastic materials (Ronchin et al. 2013). The models are implemented with a discretization of the prismatic volume ( Fig. 2b) containing the reservoir.
Due to the need of satisfying both low-velocity tomographic contrast interpreted as a shallow reservoir (Finlayson et al. 2003) and the hypocentres' location (Mori & McKee 1987;Jones & Stewart 1997;Saunders 2001), the chosen prismatic volume hosting the FEM sources at the centre of the caldera (i.e. pressurized cavities) occupies most of the low-velocity zone imaged by the tomography (Fig. 3c) and laterally encloses the hypocentres (Figs 3a and b). Keeping in mind that no shear is possible in fluids, and thus considering the hypocentres as lateral boundaries for the magma storage zone in the upper part of the reservoir, we estimate that the central prism hosts the unknown reservoir in a conservative way, although it does not entirely occupy the lateral extension of low seismicity (Fig. 3c). The prism is confined between 2000 and 5000 m depth and thus is fully included in the intrusive block. The deeper side of the prism is far below the deeper hypocentres recorded along the elliptical faults and below the bottom of the shallow reservoir imaged by the tomography (Finlayson et al. 2003), ensuring the investigation of a shallow magma reservoir constrained by the geophysical studies (Fig. 3).
A complete description of an FEM model includes the governing equations, initial conditions, boundary conditions, applied loads and the tessellated geometric domain. The 3-D problem domain approximates a solid hemisphere with a diameter of 100 km  Finlayson et al. (2003). Dashed lines on the slices at 2 and 3 km represent the seismically active portions of the ring fault above 2 km and below 3 km, respectively. Black frames represent the footprints of source array. centred on the caldera, whose upper surface is a stress-free surface represented by the topographic and bathymetric relief. The FEMs configuration is summarized in Table 2. At the outer surface of the hemisphere, far enough from the caldera for displacements to vanish, we specified zero displacements (Figs 2a and d). The governing equations describe elastic behaviour in a 3-D domain having no body forces and a spatial distribution of isotropic elastic material properties E (Young's modulus) and σ (Poisson's ratio).
The FEM domains are formed by two parts: (1) a heterogeneous part made of six geological blocks with a shape and size described in Ronchin et al. (2013) and defined by tomography and geology (Figs 2a, c and d), and (2) a central prism surrounded by the intrusive block (purple body in Figs 3a and b), both of them having the same material properties. The bulk dynamic elastic properties of the geological blocks (Table 2) are calculated from the Vp velocities provided by the tomography (Finlayson et al. 2003) and scaled to static values (see Supporting Information for further information). The equivalence between dynamic and static elastic properties is valid for depths greater than a few kilometres (Simmons & Brace 1965). The central prism is meshed with cubic elements (150 m edge length) that ensure the generation of a regular array of sources. The surrounding domain is entirely partitioned with tetrahedral elements, due to their ability to assume complex shapes (e.g. topography). On the free surface, the characteristic length of the tetrahedral element edges is about 200 m for the elements above the prism containing the sources, where a higher displacement gradient is expected, and gradually increases up to about 3000 m near the farfield boundary, where lower displacement gradients are expected. Continuity of displacements between the two different meshes of the two parts is guaranteed by coupling the displacements of element nodes along the sides of the central prism where the two parts have the same number of nodes at the same position, thus forming pairs of nodes (see Supporting Information). The central prism (Fig. 2b) is made of (1) a core from which elements are removed to generate the array of pressurized sources and (2) a shell of cubic elements not used for the source generation. Multiple FEMs, each having a different ith pressurized cavity (300 m edge dimension), are sequentially generated by removing a different group of eight elements (150 m edge dimension) from the core of the common FEM model mesh (Fig. 2b). Elements of the core not removed at the time of the ith source generation are part of non-active sources and are structural elements having material properties of the intrusive block in which they are included. Duplicate material properties characterize the elements of the shell that are not used to generate the cavities (Fig. 2b). Their role is to guarantee a smooth transition between the cubic sources and the tetrahedral mesh, and a uniform discretization of the volume around all sources. This ensures that all the sources equally behave as a centre of dilatation (i.e. McTigue finite source) (see the Supporting Information, Figs S1 and S2). All eight elements chosen to generate the ith cavity are removed only once. This means that a virtual regular array of non-overlapping cavities is formed and that we obtain 2160 cavities arranged in nine vertical layers (Figs 2b, 3a and b) from a total of 17 280 cubic elements constituting the core.
Initially, each elementary cavity is generated and filled with a compressible hydraulic fluid (i.e. magma having bulk modulus k and initial density, pressure and volume) by: (1) removing eight elements from the central prism to form an ith cubic cavity and (2) coating the cavity with Abaqus volumetric hydrostatic fluid elements (Abaqus 2009), which fill the cavity with a mass of magma in equilibrium. The source of deformation is then initiated by injecting a fluid load q (in kg) into the filled cavity through the hydrostatic fluid elements. In response to q, the current hydrostatic fluid pressure P c is given via the equation: where ρ is the reference density of the fluid at zero pressure, ρ in is the initial density of fluid filling the cavity (where ρ in = ρ, if no tem-perature dependency of the fluid is considered) and ρ c is the fluid density at current pressure. The process is isothermal and k is independent of change temperature. Properties of the fluid (i.e. magma of Rabaul reservoir) used for the computation of the Green's functions are summarized in Table 2 (see the Supporting Information  for further information). By removing the elements from a common prism, we do not have to remesh the model in order to generate new sources, thus saving time and computational resources. The final isotropic pressure change P i , as well as the volume change V i of the ith cavity of the array, are solutions variables computed in response to: (1) the fluid mass q i injected through the hydrostatic fluid elements (Abaqus 2009) and (2) the accommodation of the cavity structure. P i is thus numerically calculated by Abaqus in terms of fluid bulk modulus, cavity geometry and the material properties surrounding the cavity. The hydrostatic fluid elements, besides allowing the injection of the fluid mass, provide the coupling between the deformation of the cavity and the change in pressure, P i . We have verified that at the pressure range, cavity geometry and material properties considered, P i is linearly related to q i . By this definition, we assume that the change of pressure is entirely and only due to a concordant change of mass in the source, the flux of magma.

Linear inversion based on a family of FEMs
We generate a family of multiple FEMs whose different cavities are non-overlapping, equal volume and distributed in a regular array. This array of sources allows a relatively simple finite-difference formulation of the Laplacian smoothing operator L, to regularize the solution of the linear inverse problem. Each cavity of the array is loaded with a unit mass of magma of 1 × 10 9 kg in order to generate one pressurized source in each FEM model of the family, whose strength will be estimated during the inversion. The fact that the prism is bigger than the possible reservoir allows us to account for the boundary conditions applied to the array of sources during the inversion process and to limit the effects of the boundaries. We assume no sources outside the 3-D array and apply Dirichlettype boundary conditions (Wang & Anderson 1982;Fig. 3a). The linear relationship with the induced displacements is validated by performing an FEM-based linear inversion of Rabaul synthetic data, provided in Fig. S4 in the Supporting Information. The forward solution for the displacements at the free surface, caused by the array of expanding FEM cavities, is the superposition of the displacements generated by each of the sources. To determine the 3-D distribution of change of pressure and the flux of magma entering or leaving the sources of the array, we estimate the elements of a vector q through a weighted least-squares inversion to solve the equation: where G is a 802 by 2160 matrix of Green's functions (displacements generated by the overpressure due to the injection of a unit mass of magma, 1 × 10 9 kg, into each FEM source of the array and projected onto the LOS), d is the vector of observed LOS displacements, W is a diagonal matrix of data weights whose W i,i element is the reciprocal of the variance σ i of the ith datum d i and L is the finite-difference approximation of the Laplacian smoothing operator (Freymueller et al. 1994), whose construction is extensively described by Masterlark & Lu (2004). By weighting the data with the corresponding SD (Fig. 1d) through W, we give higher influence to the data with steady LOS displacements (assumed to be related to the overall deflation of the entire reservoir, the target of the study). The importance of L is governed by the size of a scalar damping parameter β that controls the relative importance of minimizing the roughness of the solution versus fitting the data. L is defined in a way that allows: (1) the application of Dirichlet-type boundary conditions (Wang & Anderson 1982) to q along the outer boundaries of the source array ( Fig. 3a) and (2) the assumption of no sources outside the source array. The phase ramping in the InSAR image, caused by the uncertainties of the satellite position (Massonnet & Feigl 1998), is taken into account in eq. (2) by simultaneously fitting a plane. To do so, G is appended with three column vectors xy1 corresponding to the position of the InSAR data (x and y) and a unit vector. The unknown linear parameters for the vertical projection of the plane (a, b and c) are appended to q. The value of the damping parameter β that weights the smoothing operator is determined from the trade-off analysis between the misfit of the predicted displacements with respect to the observed displacements and the roughness of the inversion solution defined by Lm 2 2 = m T L T Lm. Note that the weighted smoothing appends the weighted G matrix. The three null column vectors [000] are necessary to preserve the matrix dimensions associated with the nuisance parameters [xy1]. By performing the inversion using the value of β at the knee of the curve (Fig. 4a), we identify a weighted least-squares damped (WLSD) inversion for solutions that is a good compromise between fitting the data versus minimizing the solution roughness (i.e. complexity of the model) and satisfying the boundary conditions. In order to obtain the model solution for the distribution of mass change (i.e. the flux of the fluid magma) in the array during the time frame of data acquisition, the estimates of the sources' strength q are scaled by the unit mass of injected magma. For each ith source, the change in pressure P i , being linearly correlated to the corresponding injected fluid mass, can be calculated by scaling the FEM solution for the unit hydrostatic change in pressure P i by the correspondent source strength estimate q i : Model estimates are fully described by providing their uncertainties and the model resolution, essential to interpreting models in terms of reliability (Aster et al. 2005). In this work, the model resolution is provided through a checkerboard test. Its aim is to show how well the estimated source parameters are independently resolved and how significantly smoothed or otherwise biased our solution compares to the true (but unknown) solution.
The uncertainties associated with model estimates depend on the covariance of the data and on the way errors are mapped from the data to the model parameters. They are described by the diagonal elements of the model covariance matrix (Menke 1989): whereĜ is the first bracketed term on the left side of eq. (2).

R E S U LT S : D I S T R I B U T E D S O U RC E S A N D T H E 3 -D G E O M E T RY O F T H E R A B AU L S H A L L O W R E S E RV O I R
The strength distribution of sources in the array (Figs 5a-i) is estimated in terms of both fluid flux and related pressure changes from the reduced InSAR data through a linear least-squares inversion (eqs 2 and 3), sweeping through a range of damping values of β. The result is an L-curve of points at the knee of which the preferred solution is visually identified (Fig. 4a). The preferred damping value, β = 200, corresponds to a root mean square error (RMSE) of 0.0091 m yr −1 and is chosen in a way that the data are neither substantially overfit, nor underfit with respect to the estimated data errors. This misfit, being of the same order of magnitude as the roughly estimated accuracy of the InSAR data (0.01 m yr −1 ) used for the inversion, is acceptable. The comparison of observed versus estimated LOS displacements (Fig. 4b) shows a diagonal distribution of points mostly included in a buffer area of 1 cm, roughly representing the accuracy of the observed InSAR displacements. A uniform, broad, but small misfit is observed west of the caldera border. The highest misfits are recorded for the points inside the caldera and in particular for the areas north of Matupit Island (Fig. 4e). The solutions of the WLSD inversion (Figs 5a-i) show the distribution of fluid flux (change of fluid mass per year) in the source array. The main negative fluxes of fluid (i.e. fluid flowing out of the source) are distributed along a diagonal direction between Vulcan and Tavurvur, forming a distribution with amorphous shape (Figs 5a-i). The diagonal trend of the distribution is emphasized by a diagonal cut on its southern side (red dashed line in Figs 5e-g). Higher negative flux values are localized in two areas. A first area is located in the NE of the array (SW of Tavurvur volcanic centre), starting at 2600 m and clearly visible from 2900 m depth (Fig. 5c). The flux under this area is quite consistent over all the layers in depth (Figs 5c-i), although its negative values are higher between 3500 and 4100 m depth (Figs 5e-g) and they remarkably diminish in the deeper layer. A second area of high negative flux values, slightly elongated north-south, is located at the western side of Blanche Bay in front of Vulcan (Figs 5d-h), starting at 3200 m depth (Fig. 5d). The higher negative flux values of this area are found northeast of Vulcan at 4100 m depth, in front of Vulcan North Shore (VNS; Figs 5f and g). Distributed transitional negative flux values appear to be gradually connected at depth in these areas, starting at 3200 m (Figs 5d-h), forming a 3-D distribution with amorphous shape. The highest negative flux, −0.12 × 10 9 kg yr −1 , is reached under the first area at 3800 m depth (Fig. 5f). The total magnitude of fluid mass flowing out of the source array is −74.5 × 10 9 kg yr −1 . Values of flux vanish on the shallower and deeper layers of the source array (Figs 5a, b and i).
To investigate the model resolution, we created a synthetic distribution of unit sources in a checkerboard pattern (Figs 61-3). The checkerboard test gives an indication of how well details are resolved by the inversion. Resolution is higher at shallower levels (Figs 6a-c), but rapidly decreases at deeper levels. Squared shapes of about 1 km side length are still recognizable in the first three levels and at the sides of the fourth, while at deeper levels the resolution is lower and the elements of the synthetic source pattern are smeared. Model uncertainties (Figs 7c-g), generally reasonable with respect to the source estimates, are higher at the southeastern corner of the array.

D I S C U S S I O N
For computational reasons, complex volcano deformation sources are generally interpreted through models having simplified source geometries. The strategy proposed here allows us to integrate advanced FEMs into the inversion of InSAR data for the generation of FEM-based images of volcano reservoirs. This satisfies the need for studying and interpreting volcanic reservoirs through models with more realistic 3-D source geometries and domains that honour heterogeneous distributions of material properties as indicated by seismic tomography and geology. Results of the Rabaul case study (Fig. 5) show how the proposed strategy improves the imaging of 3-D deformation sources by identifying the location and geometry of the active parts of the reservoir. This allows tracking the evolution of reservoirs, which may thus widen the understanding of reservoir dynamics, advance our ability of monitoring active volcanoes and improve the prediction of eruptions.
The inverse solutions for a distribution of magma withdrawal (i.e. mass of fluid leaving the sources) and the related pressure change track the extension and shape of the shallow magmatic reservoir under the Rabaul caldera (Figs 5 and 8). The solutions of mass removal show a complex distribution made of two areas of high values: a first one located at the northeastern corner of the array, starting at shallower levels, merges at higher depths (between −3500 and −4400 m, in Figs 5e and f) with a second one located on the western side of the caldera, NE of Vulcan. When considering the shape of the amorphous distribution of sources as possible source of deformation, we must keep in mind the resolution of the inversion model used. In our case, the resolution of the model (Figs 6a-i) allows us to detect the presence of two areas of higher negative flux values at the borders of the source array, starting at different depths (Figs 5c and d). For depths deeper than 3500 m (Figs 5e-i), no further details than a smeared distribution of sources at the centre of the array can be achieved due to the lower model resolution (Figs 6e-i). Thus, the inverse solution favours shallow sources at the sides of the array where the resolution of the model is higher (Figs 6a-c). Higher resolution is found near and beneath the coastline, as no deformation data are available from the bay. No data coverage over the bay, together with the way the data variance is mapped to model parameters, determines higher uncertainties for source estimates in the southeastern corner (Figs 7c-g), while smaller uncertainties correspond to source estimates at the western and northern side of the array.
Correspondences between the estimated distribution of fluid fluxes and the geological features, petrographic information and observations of the volcanic activity support our solution. The elongated distribution of fluid mass withdrawal, oriented SW-NE between Vulcan and Matupit Island and centred at about 3500 m depth (Fig. 8a), does not laterally intersect the seismicity recorded at Rabaul caldera central elliptical faults between 1971 and early 1992 (Fig. 5). This is consistent with the assumption that the faults do not propagate in the molten rock and define a physical constraint for the reservoir. In particular, we observe that the distribution of fluid flux in all levels is confined to the shallower seismicity (<2 km depth, red dots in Fig. 5b), resembling its shape on the southeastern side (dashed red line in Figs 5e-h). This suggests a possible relationship between the extension of the imaged reservoir and the faults related to the shallower seismicity. The old submarine cones inside the caldera, located at the border of the imaged reservoir (Fig. 5a), suggest a possible correlation with the southernmost extension of the reservoir. Folds mapped during seismic reflection surveys (Pono 1990) show SW-NE trends on top of the imaged reservoir, resembling the elongation of the reservoir itself (Fig. 5a). The higher and shallower flux values, located SW of Tavurvur (Fig. 8b) and having broad vertical extension, are consistent with a higher withdrawal in this volume of the array due to the almost continuous vulcanian eruptions at Tavurvur. The western side of the amorphous distribution of sources is consistent with a batch of magma that could have fed the 1994 eruption at VNS (Fig. 5a). Furthermore, its elongated N-S shape, with slightly higher values on its northern part (Fig. 5f), resembles the distribution of historical cones (Vulcan Island, Vulcan and VNS) on the western side of the caldera (Fig. 5a).
By assuming that the observed deformation is caused by sources confined to the prism depth, we neglect: (1) the feeding dykes from below and (2) the very shallow processes having small deformation wavelength and short-period deformation impulses, such as localized hydrothermal processes or the intrusion of shallow dykes that connect the reservoir to the Tavurvur eruptive centre. The first point implies that we are modeling a half-opened system in which we assume that the mass of fluid can only flow out of the system, ignoring any new input of magma from deeper reservoirs. This assumption is a reasonable approximation based on the fact that the absence of petrographic changes in the volcanic products from the October 2006 eruption (RVO reports in Smithsonian Institution 2013) suggests that no significant replenishment of the shallow reservoir occurred during the period studied. The second point arises from the goal of focusing the investigation on the shape and extension of the reservoir by using its long-lasting deformation signal, rather than on shallow and short-period processes, such as the feeding of local vents or the hydrothermal activity, which is still poorly understood (Johnson et al. 2010).
The interconnection of sources at depth (Figs 5 and 8) finds correspondence in the products erupted in 2006 October and during the study period, as well as in products erupted at Rabaul during previous eruptions (Patia 2004;Bouvet de Maisonneuve et al. 2015). In fact, the andesites erupted during 2006 October are the results of the mixing and mingling of two magmas from the same system at a depth calculated between 2.7 and 3.8 km: a replenishing basalt and a resident dacite (Bouvet de Maisonneuve et al. 2015). Less evolved magmas are erupted from Tavurvur, on the eastern side of the caldera, suggesting a more direct feeding by basaltic magmas from a deeper part of the system on this side of the caldera (Roggensack et al. 1996). Dacites are the only products erupted from vents on the western side of the caldera (Roggensack et al. 1996), thus suggesting that the western area of the reservoir imaged by the solutions could represent part of the resident dacite batch. The estimated interconnection of sources at depth suggests a possible migration and interaction of fluids from opposite sides of the reservoir, which is consistent with the mixing and mingling processes observed by other authors (Patia 2004;Bouvet de Maisonneuve et al. 2015).
The existence of two sources at Rabaul has been proposed by previous analytical half-space homogenous models, based on the data collected during the seismo-deformational crisis of 1983-1985: one source located southeast of Matupit Island, at a depth between 1.2 and 2 km (Fig. 5a), and another under Vulcan at a depth of 3 km ( Fig. 5d; McKee et al. 1984McKee et al. , 1989Archbold et al. 1988). The distribution of our solutions resembles these sources, although locating them at higher depths. Thus, our study confirms and extends the results of these early studies. Compared to the previous models, the proposed model allows a more realistic domain with heterogeneous distribution of material properties and a corresponding source geometry that is determined by the deformation data. In fact, the model depicts reservoir details unachievable by using single a-priori sources and features not solved at the resolution of the tomography. Moreover, it allows detecting the connection between two relatively shallow batches of magma where the mixing and mingling of the two magmas can occur, as already hypothesized by petrographic studies (Patia 2004;Bouvet de Maisonneuve et al. 2015).
The estimated total amount of mass withdrawal per year is 74.5 × 10 9 kg yr −1 . This corresponds to an estimated total volume of 119 × 10 6 m 3 , withdrawn by almost continuous Tavurvur vulcanian eruptions during the entire study period. Although we do not have independent information about the erupted volume to validate our estimate, we believe that we are proposing a reasonable value, considering that the volume of the 2006 eruption, 200 × 10 6 m 3 (see Saunders et al. 2007), could be a reference value for a typical eruptive volume over a short period at Rabaul.
Pressure values calculated for the sources (Fig. 5) are realistic, being smaller than the shear strength of diorite in which the sources are located. This, at first approximation, prevents a slip along fractures. In fact, no collapses that were observed during the period studied.
A little area of anomalous positive LOS velocities is observed along the southwestern side of the old airport strip, north of Matupit Island (Fig. 4e). This anomaly may reflect a local groundwater or soil compaction signal, since this is a drainage area. The small but broad misfit west of the caldera could be generated by tectonic movements or deeper sources not considered in our model.
The distribution of pressure changes, provided as variable solutions by the hydrostatic elements (Abaqus 2009), is valuable information for studies whose aim is the understanding of the magmatic processes occurring inside a reservoir, which generate pressure changes responsible for surface deformations. In this work, a principal assumption is that the change of magma mass, which changes both the pressure and volume of a given source, is entirely due to a magma mass flux into or out of the system. It is important to note, however, that alternative mechanisms could be responsible for the pressure and volume variation (e.g. crystallization, degassing, remelting, expansion of the hydrothermal system, etc.) that may have different relationship with the change of mass in the system. In further studies, with gravity data available, the gravity signal could be modeled and a fully joint inversion of gravity and deformation could be performed in order to discriminate the processes responsible for the deformation. Therefore, we believe that the proposed model, by simulating change of mass, pressure, volume and deformation, could open up further studies combining deformation and microgravity data that could better shed light on the mechanisms responsible for surface displacements at active volcanoes.

C O N C L U S I O N S
In this work, we propose a computationally inexpensive strategy for an FEM-based linear inversion of InSAR data capable of estimating the extended, arbitrary distribution of mass and pressure changes within the integration domain of FEM models. This method leads to the creation of 3-D images of the sources responsible for surface deformation. The method pushes FEM-based inverse modeling a step forward by removing the requirement for a-priori assumptions regarding the shape of volcano deformation sources and letting the data dictate the shape of the deformation source as a spatial distribution of mass flux, and thus a distribution of pressurization throughout a geometrically complicated reservoir.
By applying the method to InSAR data of Rabaul caldera, we have been able to image the extent and shape of the active reservoir (Fig. 8b). Our study images the presence of two main sources at opposite sides of the caldera and connected at depth, showing that the shallow magmatic system is an interconnected reservoir instead of distinct magma bodies. This implies a high mobility of magma between different parts of the system and explains how a long-lasting eruption at one side of the caldera, in this case at Tavuvur, can trigger a long-term deflation of the entire system.

A C K N O W L E D G E M E N T S
ER thanks Tobias Meinel, Carlo Fonda, Gabriella Bacchelli and Sergio Sokota for their constructive feedbacks and their help in proofreading the manuscript. ER thanks Maria Charco and Joachim Gottsmann for their valuable suggestions and Caroline Bouvet de Maisonneuve for elucidations about the petrography of Rabaul volcanic products. ALOS data were provided by JAXA with access provided by Geoscience Australia. TM's contribution was supported, in part, by NSF-EAR 1316082, 0943943 and 1264290, and NASA Earth Surface and Interior NNX17AD96G. The study was conceived by ER and designed by ER and TM; ER performed the research; JD processed the InSAR data; ER and TM analysed the results; SS and JM provided help for the interpretation of the results; ER, TM and JD wrote the manuscript. We thank the Editor, Prof Duncan Agnew, Maurizio Battaglia and an anonymous reviewer for their constructive comments that greatly improved this paper.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJIRAS online. . Details of cavity discretization of each tested source are shown: cavity types L1 and S1 made by removing one element, cavity type L8 and L8i made by removing eight elements (colours as in panel b for location reference).    Table S1. Elastic formulae, size and elastic properties of the 3-D geological bodies. Table S2. Values of parameters used to calculate the bulk modulus of Rabaul magma.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.