Determination of Cerro Blanco volcano deformation field using method of fundamental solutions

ABSTRACT According to geodetic research works, surface deformation in the form of uplift or subsidence in volcanic areas is either a sign of magma moving towards the opening of the volcano (inflation) or removal of the magma source (deflation). Using the new method of fundamental solutions (MFS) in this study, a deformation field for the surface of the volcano is calculated considering the effect of topography. MFS is a numerical technique for solving boundary value problems with known partial differential equations. This technique has not been used in volcanic deformation studies so far. Because of the simplicity and efficiency of the technique, it is also an effective tool for solving a wide range of problems in other fields of science and technology. To test the method, the displacement calculated using the MFS was compared with that of the interferometric synthetic aperture radar observations from the previous study in Cerro Blanco volcano. The volcano was in the deflation mode during this period at the rate of 1.2 cm/yr. The comparison showed a root-mean-square error (RMSE) in the order of 2 mm which represents a satisfactory agreement with the results of the observations, less than the RMSE of the analytical models considered.


Introduction
Volcanic eruptions often cause irreparable damages to human life. During eruption, high-intensity magma rises towards the surface. Geodetic observations indicate volcanic uplift or subsidence in many parts of the earth which is attributed to the transfer of magma at depth (Mogi 1958;Davis 1986;Dvorak & Dzurisin 1997). Understanding the geometry of the magma reservoir and ground deformation in volcanic regions is crucial in predicting volcanic eruptions and prevention of adverse events. Many studies have been conducted in volcanic and high-risk geodynamic zones like Yazdanparast & Vosooghi (2014), Jebur et al. (2015), and Deffontaines et al. (2016). One of the pioneering studies in this area was done by Omori in 1910 who calculated the vertical deformation of a volcano in Japan using the precise levelling (Omori 1914). Mogi calculated surface deformation of a volcano using the parameters for its magma reservoir such as position and depth (Mogi 1958). Surface deformation was measured before and after the eruption of Mount Etna in 1981 using precise levelling and gravimetry (Carbone et al. 2008). Electronic distance metre (EDM) technique was also used to calculate the deformation in Etna (Carbone et al. 2008). With the recent progress in the application of global positioning system (GPS) and interferometric synthetic aperture radar (InSAR) in geodetic observations, many studies have been conducted in this area. Development of volcano deformation models has been revolutionized over the past two decades. Davis (1986) developed Mogi Model for a non-isotropic point source. Yang et al. (1988) presented a solution to a hole being uniformly formed under pressure in an elastic half space.
These solutions have been successfully applied to interpret the deformation of resources in Wu and Wang (1988), and Fialko et al. (2001) investigated the deformation due to a flat circular source in an elastic half space. Beauducel et al. (2004) presented a model for Campi Felegri volcano using the boundary element method (BEM). Di Traglia et al. (2015), Singhroy et al. (2004), Shirzaei et al. (2011) and Brunori et al. (2013) studied the deformation of Stromboli, Miayke-Jima, Damavand and Cerro Blanco volcanic regions, respectively, using InSAR observations. Masterlark et al. (2010) developed finite element models (FEMs) of Okmok volcano to simulate its deformation. Charco and del Sastre (2014) presented a procedure for FEMs which can be combined with explorative inversion schemes in volcanoes.
The effect of topography in two dimensions has been considered by Mctigue and Stein (1984) and McTigue and Segall (1988). Cayol and Cornet (1998) have also investigated the effect of topography on deformation models at Mount Etna using a BEM model. Lungarini et al. (2005) used a numerical finite element technique to model rigorously the three-dimensional (3D) topographic effect of ground displacement at Mount Etna.
Using the new Method of Fundamental Solution (MFS) in this study, a deformation field for the surface of the volcano is calculated considering the effect of topography. MFS is a meshless boundary method applicable to certain boundary value and initial/boundary value problems. In particular, it is readily applicable to homogeneous boundary value problems in which a fundamental solution of the operator in the governing equation is known explicitly (Karageorghis et al. 2011). MFS, since its introduction as a numerical method by Mathon and Johnston (1977), has become increasingly popular and been used for solving various types of problems by Marin and Lesnic (2004), Johansson and Lesnic (2008), Ahmadabadi et al. (2013), Reeve (2013), Ala et al. (2015), Wang and Zheng (2016), and Shidfar and Darooghehgimofrad (2017).
The popularity of the MFS is primarily due to the ease of its application in various geometries including 3D problems. The advantage of the MFS over the more classical domain or boundary discretization methods and the reasons why its application is straightforward are explained herein. MFS is meshless and a collection of points is only required for the discretization of the problem. It does not involve integration which is potentially troublesome and complicated, as is the case with the BEM, for example. Finally, it is a boundary method which shares all the advantages that the BEM has over domain discretization methods such as the finite element (FE) or finite difference (FD) methods. Therefore, MFS is ideally suited for the solution of problems in which the boundary is of major importance or requires special analysis (Karageorghis et al. 2011). Finally, in mesh-based numerical methods, human intervention can never be completely eliminated from the mesh generation process (Mossaiby et al. 2016). This is evidently not the case in the meshless MFS method.
An introduction to the MFS and its application in deformation field of volcano is presented in the following section. The third section includes a brief presentation of Cerro Blanco volcano as the study area and the associated data-set used in this study. The fourth section contains simulation tests. The results of implementation of the MFS on Cerro Blanco will be presented in the fifth section and, finally, the conclusion of the investigation is discussed in the last section.

Method of fundamental solutions
MFS is a method for the solution of certain elliptic boundary value problems. The solution in this method is approximated by a set of fundamental solutions of the governing equations expressed in terms of sources located outside the domain of the problem. Besides the advantages of the MFS mentioned in the previous section, power sources are located outside the domain of the problem, so the solution will not encounter the obstacle of singularities which is the case in the BEM. The MFS is a numerical method for solving partial differential equations (PDEs) when the fundamental solution of the PDE is known. To illustrate how the MFS is implemented, assume a bounded domain V with boundary G as shown in Figure 1 and a PDE denoted by the following equation: where L is a linear differential operator with constant coefficients and u is the unknown parameter which should be estimated, such as displacement, stress or a combination of both.
To solve these equations, one of the boundary conditions in the form of Dirichlet, Neumann or Robin should be considered.
An approximation which satisfies both the governing PDE and the boundary condition (at a finite number of points) is constructed in the MFS. The fundamental solution G(P, Q) of Equation (1) satisfies the following equation: where d is the Dirac delta function, Q is called a source point and P is the boundary point. As shown in Equation (2), fundamental solution is the system response when the simulator is the Dirac delta function.
If the source point Q is placed outside the domain V, then the fundamental solution G(P, Q) and any linear combination of fundamental solutions satisfies Equation (1) without any singularity. Hence, an approximation to Equation (1) is constructed in the form of where Q j represent N source points placed outside the domain V and c j are constant coefficients. u M is the boundary value at point P. Figure 1 is a schematic representation of the boundary and the source point's distribution. To determine constants c j , that it is performed with the use of boundary points (blue circle in Figure 1) and taking into account boundary conditions for Equation (3). Assuming M boundary points and N source points, there is a system of equations with M equations and N unknowns. Now that the constants have been determined, the value of u at any point within the domain can be calculated using suitable fundamental solution between that specific point and the source points.

Application of MFS in deformation field of volcano
Surface displacement in and around topographic configuration comes from the interaction of the ground surface and perturbed displacement field (Luz on et al. 1997). Hence, the total motion may be explained as the superposition of the free field u f i , which corresponds to the solution for volcanic loads (pressure changes) in the half space, and the field perturbed by the ground topography, u t i : u f i is usually computed from analytical relations depending on the type of magma source (point, ellipsoidal, etc.), hence the MFS is used to determine the perturbed field. In the absence of body forces, the governing equations of equilibrium for a homogeneous isotropic linear elastic solid are the Cauchy-Navier equations. Using the indicial tensor notation in terms of the displacements u t i , i = 1, 2, 3, the Cauchy-Navier equations, in a bounded 3D domain V of the solid, take the dimensionless form of (Poullikkas et al. 2002 where λ and m are the Lam elastic constants. Summation over repeated subscripts is implied and partial derivatives are denoted by Boundary conditions of Equation (5) are considered in the form of The operators B i specify Dirichlet, Neumann or Robin boundary conditions and t i , i = 1, 2, 3, are tractions. The perturbed displacements are approximated by linear combinations of fundamental solutions: If available boundary values are tractions, their fundamental solution can be calculated from the fundamental solution of displacement. In the above equations, N is the specified number of sources, G ij , i, j = 1, 2, 3, are fundamental solutions that can be computed for a full-space from the following equation: where R 2 = (q 1 ¡ p 1 ) 2 + (q 2 ¡ p 2 ) 2 + (q 3 ¡ p 3 ) 2 , R i the partial derivative R with respect to boundary point p i , n is the Poisson ratio and p i and q i , i = 1, 2, 3, are the coordinates of boundary and source points, respectively. Fundamental solutions for the half space are obtained from (10) (Mindlin 1936).
Because of the symmetry of the problem, the displacements caused by force acting along the second axis are calculated by replacing q 1 ¡ p 1 and q 2 ¡ p 2 components in G i1 : (8) are vectors containing unknown coefficients, Q contains the coordinates of the sources Q j , which lie outside domain and a set of points of p i , i = 1, …, M are boundary points (Poullikkas et al. (2002)).

Study area and used data
The study area is Cerro Blanco volcanic field between 26 39'S-26 54'S latitudes and 67 36'W-67 54'W longitudes in the Catamarca Province in Argentina ( Figure 2). It is a volcano collapse structure located at an altitude of 4400 m. The caldera has been active for the last eight million years and eruptions have caused several ignimbrites. One of the most recent eruptions occurred 73 ka. About 5 ka, the largest volcanic eruption of the Central Andes emerged at Cerro Blanco, creating the most recent caldera as well as thick ignimbrite layers. Cerro Blanco has a high desert climate with high insolation and long-term aridity. Plant growth is minor, so it is suitable for the application of InSAR technique. Principal characteristic of the Cerro Blanco is its uniqueness in terms of deformation among the actively deforming volcanoes as it is subsiding (Pritchard & Simons 2004;Henderson & Pritchard 2013). Limited research has been carried out specifically on deformation and structure of this volcano (e.g. Brunori et al. 2013;B aez et al. 2015). To fill in this gap, the deformation condition of this volcano is investigated and presented in this study. Weather conditions and absence of vegetation are beneficial in using InSAR observation in this volcano. The data used here are multiple radar images from the European Space Agency satellite ENVISAT. The data-set includes 14 images in descending orbit (track 10), and spans the period 2003-2007. Images are in VV polarization, with an incidence angle of 23 and azimute of -168 . The processing chain in StaMPS-MTI method previously was implemented by Hooper (2008). StaMPS-MTI allows exploiting two A-InSAR (Advanced InSAR) techniques, the PS technique (Ferretti et al. 2001) and the SBAS technique (Berardino et al. 2002). Combining these two approaches, the deformation signal for higher number of points is acquired with an overall signalto-noise ratio higher than that of each single technique separately. This leads to more reliable results, in particular, in cases with limited data-set such as the one used in this study. Mean velocity fields and the deformation history of the Cerro Blanco were then calculated. The retrieved map shows a distribution of negative and quasi-circular deformation pattern (movement away from the satellite in the line of sight [LOS]). The measured rate of subsidence has a maximum value of ¡1.2 cm/yr in the radar LOS (Brunori et al. 2013).

Simulated tests
Before using the method on a real volcano, it was tested for validation in simulation mode without considering the topography effect. To organize the simulated data-sets as the boundary values, spherical (Mogi 1958) and ellipsoidal (Yang et al. 1988) sources were considered. In reality, observations are tainted by random errors. Therefore, in the current test to consider the effect of random noise, a normally distributed component with a variance of 0.5 cm 2 to the boundary values is added. This variance is comparable to the real variance of modern radar interferometric methods (Berardino et al. 2002). For Mogi source type, a square as the surface of half space with dimensions 20 km £ 20 km with a spherical cavity located inside it as the inner border was modelled simultaneously. The parameters chosen for the cavity are shown in Table 1. The MFS was applied by taking 25 boundary points on the outer boundary as well as 25 points on the surface of the spherical cavity and the same number of source points in parallel to the boundary points outside the outer boundary and inside the inner boundary. These points were evenly distributed as shown schematically in Figure 3. For ellipsoidal type, an ellipsoidal cavity instead of a spherical one was considered as the inner border.
The corresponding parameters are presented in Table 1. Source and boundary points are similar to spherical cavity. Geometry of this cavity type is illustrated in Figure 4. To obtain the amount in the outer boundary points, the relationship introduced by Okada (1992) for an inflation point source and Yang et al. (1988) for ellipsoidal cavity were used.
The traction boundary values for the cavity similar to Karageorghis et al. (2011) were assumed zero. Poisson ratio and shear modulus were assumed as 0.25 and 7 GPa, respectively. Fundamental solutions for half space are obtained from Equation (10).
After calculating the fundamental solutions to displacement and tractions, the coefficients a, b and c are captured as described in Section 2.1. Because one of the typical observations for detecting the Earth deformation is InSAR data, simulated movements of the outer boundary were projected along the LOS of satellite descending orbit using the formula was presented in Hanssen (2001) similar to the radar data. Displacements were sampled in uniform intervals 0.1 km, until these are comparable with the spatial resolution of InSAR displacements.
To interpret the movements derived from the MFS, the displacements of the spherical (Mogi (1958)) and ellipsoidal (Yang et al. (1988)) volcanic reservoir models were applied. After the projection Table 1 of these movements along the LOS of satellite, the differences were calculated as shown in Figures 5  and 6. The difference amounts show that a satisfactory agreement between the modelled displacements and displacements resulted from the fundamental solution. The differences are within §3 mm. It should be noted that before the selection of the input parameters in the MFS shown in Figures 5  and 6, the influence of two important parameters of MFS were investigated. These include the distance of source points to domain and the number of boundary points that the effect of its changes was studied.
To check the influence of the former parameter (i.e. the distance of source points from the domain), the problem was solved for several different modes while the other parameters were fixed.
The results are provided in Figure 7. As shown, the error value decreases with the increasing distance. Because of the distance of 5 km to the next, there was no significant change in the amount of error, this is intended for points. The second parameter was the number of boundary points that the effects of its alterations were studied. Again, other parameters were kept fixed while investigating the effect of number of boundary points. The results are shown in Figure 7. It is evident that with increasing the number of the boundary points, the error decreases.
Due to the fact that with boundary points over 25 on each border, there was no significant change in the amount of error; this number of points is selected. In order to investigate the effect of topography on deformation, a numerical test was carried out by Williams and Wadge (1998) using 3D MFS. They demonstrated that topography has a significant effect on predicted displacement, particularly for magma chambers at relatively shallow depths (depths in the order of topography height). On this basis, a source located at shallow depth compared to the topography was considered in this study. Figure 8 represents the geometry of the model.
The effect of topography is illustrated by volcanoes of altitude H = 0, 881, 1819 and 2886 m and average slope of the flanks of a = 0 , 10 , 20 and 30 . The volcano edifice is a cone with the radius of 5 km. Poisson ratio and shear modulus and magma chamber parameters are similar to the previous tests. Hundred boundary points were distributed on the cone uniformly. Displacements perturbed by topography on the boundary points are obtained from the free-field displacement, u f i , caused by a centre of dilatation and the method of Williams and Wadge (1998) assuming a different magma source depth at each point for which a solution is required. For fundamental solution of displacement in the topography space, Equation (9) was used. Figure 9 represents displacement changes caused by a spherical source located in an elastic homogeneous medium. The results of changing surface slope, a, and height variations, are compared with flat free-surface solution, a = 0 , H = 0 m. Topography has a significant effect on the magnitude of displacements. With increasing height especially near the symmetry axis, vertical displacements decrease in magnitude. For the surface slope of 30 , the summit is a place of local minimum while the maximum radial displacement emerges at about 2 km from the centre of dilatation. In radial displacements, these reductions near the summit are smaller. Displacements in LOS direction have almost similar trend to vertical  displacements. Cayol and Cornet (1998) pointed out that the interpretation of ground surface displacements without considering topography can lead to erroneous estimations. They found that the steeper the volcano, the flatter the vertical displacement field. The curves corresponding to vertical displacements are slightly flatter in the proximity of the symmetry axis. The location of maximum value for radial displacement varies depending on the height of summit. Vertical displacement caused by a centre of dilatation diminishes faster at large horizontal distances (r > c) when topography is neglected.
These results are in a satisfactory agreement with those of Cayol and Cornet (1998). The difference is that in this study the meshless method of MFS was used instead of BEM. The advantages of this meshless method over similar numerical methods were mentioned previously.

Application of the MFS on Cerro Blanco volcano
In addition to the advantages mentioned in the previous sections, 3D MFS unlike methods such as FEM does not needs a 3D mesh generation which is very time consuming. A 3D MFS is, therefore, generated in the case of Cerro Blanco volcano in this study. A digital elevation model (DEM) provided by shuttle radar topography mission (SRTM) as topography is used in this study (Figure 10).
The implementation method is similar to that of a simulation test. A system of Cartesian coordinates with the origin located at the sea level, just below Cerro Blanco summit, is considered. In this   case, xand y-axes are orientated along WE and SN directions, respectively, whilst the z-axis points up out the medium. To select the magma source parameters, genetic algorithm (GA) was employed. Upper and lower bounds in the search space for magma chamber centre position is equal to InSAR data coverage, source depth is between 4 and 12 km within the ranges of suggested by Pritchard and Simons (2004) and Henderson and Pritchard (2013). Equation (12) allows to define a lower bound for volume change, Dv, which can be obtained from the relationship between volume change at the source in elastic half space and LOS displacement (u LOS ) observed for a Mogi source accordingly : which is À9 Â 10 À5 km 3 . The Mogi source volume changes lie between ¡9 £ 10 ¡5 km 3 and À9 Â 10 À3 km 3 . Table 2 shows the parameters used in this modelling. The volume change selected is about 20% less than that estimated by Henderson and Pritchard (2013). The latter is based on the average local topography. Four hundred points were selected on the boundary to form a uniform grid. Poisson ratio like other volcanoes was assumed as 0.25 and the shear modulus similar to Ruch et al. (2009), at 20 GPa. According to Cayol and Cornet (1998), the depth of magma reservoir used for prominent volcanoes, must be considered from the summit. So, the depth of magma source was assumed 5.6 km below the sea level for Cerro Blanco in this study. In order to check the results, after the implementation of the MFS, displacement field with spatial resolution of 0.1 km similar to the displacement of radar observations was collected and projected along the LOS orbit.Observational displacement field and displacement field of the MFS are shown in Figure 11(a,b) respectively. As shown in these figures, the results of the MFS are consistent with the observed displacements showing a subsidence rate of 1.2 cm/yr. RMSE is small and about 2 mm. In this case, subsidence has a limited elongation in the WE direction and is decreasing with the distance from the crater. Most movements occurred in the vicinity of the summit. Usually volcano displacement field is calculated using the Mogi analytical model and without considering the effect of topography. In addition, a common procedure for considering the effect of topography is to choose a constant reference height above the sea level. To compare the values obtained from these methods with the method presented in this study, three separate modes, namely without considering the effect of topography, the reference height of 2200 m (somewhere between the sea level and summit) and 4400 m (equal to summit) were considered. The displacements corresponding to each mode are shown in Figure 11(c-e), respectively. Source and elastic parameters are similar to those of the MFS method. All deformations were projected along the LOS similar to InSAR observations. Figure 11(c) shows displacement field in which the effect of topography is not considered. As shown, movements in this case are three times larger than the displacements resulted from the MFS. These movements are concentrated in a small area around the summit.This is because displacement is a function of magma source depth and that the source depth should be considered from summit, so it is located 5600 m below the sea level. Overestimation of deformation has occurred in this case and the RMSE is 6.34 mm. Figure 11 (d) displays the results in case that the reference elevation is assumed 2200 m above the sea level. This model adds a constant elevation as the average topography to the free surface. Displacements are less overestimated than the last mode (no topography). This overestimation is about two times larger than the MFS results. RMSE in this case is 3.63 mm. In Figure 11(e), movements are shown when the reference elevation is considered equal to the Cerro Blanco summit. The results in this case are underestimated with respect to the MFS displacements. RMSE in this case is 2.53 mm. Based on the cases investigated in Figure 11, it can be concluded that the effect of real topography in determining the Figure 11. (a) InSAR LOS displacement (Brunori et al. 2013). (b) MFS displacement resulted from a magma source located at 10 km depth below the summit or 5.6 km under the sea level. (c) Mogi displacement field when free surface is at the sea level and topography completely neglected. (d) Analytical displacements resulted from a magma chamber within a half space when reference elevation at 2200 m above the sea level (somewhere between the sea level and summit) is considered. (e) Analytical displacements resulted from a magma chamber within a half space when reference elevation at 4400 m above the sea level (equal to Cerro Blanco height) is considered. displacement field for high volcanoes such as Cerro Blanco is remarkable. Since the application of the MFS method is computationally simple and not very time consuming, this method can easily be used to calculate the deformation field of volcanoes. In order to analyse the deformation caused by the volcano, strain tensor was calculated. Since the components of the strain tensor are dependent on the coordinate system, invariant scalar quantities such as maximum shear of the tensor are extracted for physical interpretation, as shown in Figure 12. According to Dzurisin (2007), components of E 13 ,E 23 ,E 31 and E 32 of strain tensor for free surface equal zero. Maximum values of strain tensor components are in the range of 10 ¡7 . The results of maximum shear are in the order of 10 ¡6 . Shear strain is a key deformation measure in understanding physical processes of the earth surface such as volcano and earthquake (Voosoghi 2000). As expected, the highest values of maximum shear can be seen near the summit of the volcano. Its maximum is about 3.4 £10 ¡6 /yr.

Conclusions
Considering the effect of topography in this study, the MFS was used to calculate the displacement field of the earth in a volcanic area. The study shows that the MFS is relatively easy to implement, does not require a lot of data preparation and, unlike BEM in which the singularities or source points are outside the domain, it does not face with the troublesome integrations. Before the implementation of the method on a real volcano, the method was tested on synthetic volcanoes with spherical and ellipsoidal magma chambers and the results were compared with analytical models. The results showed differences in the range of 3 mm which is in a satisfactory agreement with those of the analytical models. The effect of topography was then shown on a deformation field. It was shown that in elevated volcanoes neglecting the effect of topography leads to large errors. Cerro Blanco volcano was selected to implement the MFS on a real case because it is an elevated volcano which can cause irreparable damage to human lives in eruption in which limited research has been done specifically on deformation and structure of this volcano. Thus, the MFS method was applied in this study to this volcano with a magma spherical source assumed to be located 10 km below the summit. The displacement field results calculated from this modelling exercise was then compared with the displacements from radar observations. Results with the RMSE of about 2 mm suggest that MFS can be relied upon as an efficient method used in investigating ground deformation. It was shown that the geometry and the deformation estimated for the volcano were different from those estimated using the analytical solutions. The movements decrease because magma source depth increases with the elevation of volcano. The deformation for three scenarios of no topography, a reference elevation somewhere between the sea level and the summit and a reference elevation at the volcano summit were calculated in this study. These led to the overestimation or underestimation of deformation with higher RMSE than those calculated from application of the MFS. Considering the effect of real topography and applying the MFS method gave realistic results. Also, with the calculation of maximum shear in the study area, it was proved that maximum shear occurred in close proximity of the summit. Cerro Blanco volcano was, therefore, proved to be a major source of surface deformation in the region. The study summarized in this paper is a preliminary research in the application of MFS in geodynamic field. Further case studies in fault zones using the real topography data and employing the MFS are recommended. Various real-time data such as those obtained from tiltmeter and geochemistry will also enhance interpretation of the surface deformation of volcanoes. Concentrating on the MFS for problems with non-Cauchy-Navier differential equation is also recommended in future studies.

Geolocation information
The study area is Cerro Blanco volcanic field between 26 39 0 S-26 54 0 S latitudes and 67 36 0 W-67 54 0 W longitudes. Cerro Blanco is a caldera in the Catamarca Province in Argentina. Part of the Central Volcanic Zone of the Andes, it is a volcano collapse structure located the altitude of 4400 m. The Central Volcanic Zone of the Andes is an area between 14 and 28 southern latitude.