Data Assimilation in Volcano Deformation Using Fast Finite Element Analysis with High Fidelity Model

Abstract

Estimation of the inner state of volcanoes are important for understanding the mechanism of eruption and reduction of disaster risk. With the improvement in observation networks, data assimilation of internal magma state using time-history crustal deformation data observed at the surface is expected to be suitable for solving such problems. Using finite-element methods capable of modeling complex geometry is desirable for modeling the three-dimensional heterogeneous crust structure, and nonlinear time-history analysis is required for considering the change in material properties due to the movement of magma. Thus, many cases of large-scale finite-element analysis is required, and the computational cost incurred is expected to become a bottleneck. As a basic study towards data assimilation of internal magma state considering change in material properties of the crust, we demonstrated that many case analyses of volcano deformation problems can be conducted in a reasonable time frame by development of a crustal deformation analysis method accelerated by GPUs. For verification of the data assimilation method, we estimated the magma trend in an actual three-dimensional heterogeneous crust structure without temporal change in material properties. We confirmed that the magma movement trend can be reproduced using the model considering crust heterogeneity, while models disregarding three-dimensional crust structure resulted in wrong estimations. Thus, we can see that using finite-element methods capable of modeling three-dimensional heterogeneity for crustal deformation analysis is important for accurate magma state estimation.

1 Introduction

Estimating the state of magma inside volcanoes are important for understanding the physical mechanism of volcanoes and predicting eruptions for disaster mitigation purposes. Available observation data is increasing with the improvement of observation networks by advance in sensor technology. For example, ground surface deformation data is available by Global Navigation Satellite Systems (GNSS) and Interferometric Synthetic Aperture Radar (InSAR) observations conducted by a variety of satellites. GNSS data have high temporal resolution, while InSAR data have high spatial resolution; thus, we can expect obtaining high resolution crustal deformation information by using both of these observation data. On the other hand, there is room for improvement in the analysis methods used for inner state estimation using these observation data. Using time-history observation data of crustal deformation, we target to estimate three-dimensional movement of magma in the crust. Here, the crust has three-dimensional heterogeneity with material property changing with the movement of magma. For such problems, data assimilation of magma state based on crustal deformation analysis using finite-element analysis which is capable of considering three-dimensional heterogeneity is expected to be effective. However, the finite-element analysis cost becomes huge and thus realization of such analysis is considered challenging. In such problems, cost for generating finite-element models with $10^{6-8}$ degrees of freedom and analyzing its linear response for $10^{3-5}$ cases is required. Thus, most methods approximate the heterogeneous three-dimensional crust with a homogeneous half space in crustal deformation analysis of volcanoes [10]. Crustal deformation analysis based on the finite-element method is started to be used; however, the scale of the analysis remains up to $2\times 10^{5}$ degrees-of-freedom for 20 cases [13]. We can see that many-case finite-element analyses required for reflecting the change in material properties is challenging.

Fast finite-element solvers capable of solving island-scale crustal deformation analysis with fault slip at inter-plate boundaries have been developed [3, 14]. We may be able to reduce the computational cost for volcanic crustal deformation by extending this method. Thus, we improve the large-scale seismic finite-element solver for application to volcano problems, and verify its applicability in this study. Together, we measure the computational cost required in many-case analyses. Furthermore, using the developed method, we conduct data assimilation for a hypothetical problem with a three-dimensional heterogeneous crust structure of an actual site to verify the estimation accuracy of time-history underground magma distribution estimation. This paper targets basic study on development of a method to estimate the movement of magma in the crust by using time-history crustal deformation observation at the surface, with appropriate modeling of the time-history change of the three-dimensional heterogeneous crust structure. The main target is development of a crustal deformation analysis method for volcanoes, and confirmation of its validity for many case analyses and applicability to data-assimilation. The change in material property of crust is slow in the target problem; thus, the characteristics of the problem is known to be constant in terms of the data assimilation. Considering the above, we neglect temporal change in crust material properties in the hypothetical problem in this study.

The rest of this paper is organized as follows. In Sect. 2, we explain the data assimilation method used in this study, and then explain the fast finite-element method developed for volcano problems. In Sect. 3, we verify the accuracy of the solver and measure the computational cost. In Sect. 4, we validate the method using a hypothetical problem with three-dimensional heterogeneous crust structure of an actual volcano Shinmoedake using GNSS and InSAR data. In Sect. 5, we summarize the paper with discussions on future prospects.

2 Methodology

2.1 Data Assimilation

We model the change in excessive pressure and the resulting expansion and contraction of the magma chamber using a spherical pressure source. In previous studies [12], it has been shown that the surface response for a spherical source can be approximated using a point source if the sphere is located in a depth more than three times its radius. In a point source, magma can be expressed as body forces. Thus, if the magma source is small compared to its located depth, we can express the distribution of excessive pressure using a body force density distribution without considering the material properties of the magma. Here, the excessive pressure of the magma source can be expressed as a linear combination of basis spatial distribution functions $\varvec{B}_k(\varvec{x})$ as

$$\begin{aligned} f(\varvec{x},t) = \sum _{k=1}^M \varvec{B}_k(\varvec{x})c_k(t). \end{aligned}$$

(1)

Here, $\varvec{x}$, t denote spatial coordinate and time, M is the number of basis functions, and ${c}_k(t)$ is the coefficient for each basis function. Using Eq. (1), we model the surface displacement $\varvec{u}_r(\varvec{x},t)$ as

$$\begin{aligned} \varvec{u}_r(\varvec{x},t) = \varvec{G}_r(f(\varvec{x},t)) + \varvec{e}, \end{aligned}$$

(2)

where $\varvec{G}_r$ is the Green’s function of surface displacement for excessive pressure at magma source, and $\varvec{e}$ is the observation error. Here, r denotes the coordinate axis, where $r=1,2,3$ corresponds to the North, East, and vertical directions. By neglecting abrupt volcanic activities such as eruptions or change in material properties due to movement of magma, we only need to consider long term quasi-static deformation for analyzing crustal deformation due to underground magma activity. Following this assumption, we approximate the crust as a linear elastic body. Due to the linear properties of Green’s functions, applying Eq. (1) to Eq. (2) leads to

$$\begin{aligned} \varvec{u}_r(\varvec{x},t) = \sum _{k=1}^M \varvec{H}_{rk}(\varvec{x})c_k(t) + \varvec{e}, \end{aligned}$$

(3)

where

$$\begin{aligned} \varvec{H}_{rk}(\varvec{x})=\varvec{G}_r(\varvec{B}_k(\varvec{x})). \end{aligned}$$

(4)

In this study, we assimilate data observed by GNSS and InSAR. Three component displacement is obtained at each GNSS observation point; thus, the observation vector $\varvec{y}_1$ obtained by GNSS is expressed as

(5)

where $t_j$ is the time at j-th time step, and $\varvec{x}_{n_1}(n_1=1,...,N_{1}$) are the coordinates of the GNSS observation points. This observation vector consists of $3\times N_1$ displacement components. We can rewrite Eq. (4) by substituting $\varvec{x}=\varvec{x}_{n_1}$ as

$$\begin{aligned} \varvec{H}^1_{ik}=\varvec{G}_r(\varvec{B}_k(\varvec{x}_{n_1})). \end{aligned}$$

(6)

On the other hand, the InSAR data only consists of the difference between two observations, with only a single displacement component parallel to the view direction. Here, we assume that we obtain one displacement component at each observation point via preprocessing. Thus, the observation vector $\varvec{y}_2$ obtained by InSAR is expressed as

(7)

where $\varvec{x}_{n_2} (n_2=1,...,N_{2}$) are the coordinates of the InSAR observation points. This observation vector consists of $N_2$ displacement components. Here, $\varvec{n}_r^T=(\cos \theta \sin \phi ,\sin \theta \sin \phi ,\cos \phi )$ is the normal vector parallel to the observation direction, where $\theta$ is the angle from the vertical direction to the satellite observation direction, and $\phi$ is the counter clockwise angle from the moving direction of the satellite measured from the North direction. In the case of InSAR, Eq. (6) is expressed as

$$\begin{aligned} \varvec{H}^2_{ik}=\varvec{n}_r^T\varvec{G}_r(\varvec{B}_k(\varvec{x}_{n_2})). \end{aligned}$$

(8)

By defining a combined observation vector $\varvec{y}=\{\varvec{y_1},\varvec{y_2}\}$,

$$\begin{aligned} \varvec{y}_i(t_j)=\sum _{k=1}^{M} \varvec{H}_{ik} c_k(t)+ \varvec{e}. \end{aligned}$$

(9)

where we define a combined observation matrix $\varvec{H}=\{\varvec{H}^1, \varvec{H}^2\}$. Although GNSS and InSAR observation data includes observation errors (e.g., local GNSS benchmark motion, GNSS reference frame errors, InSAR planar correction) that must be considered during data assimilation, we assume these are removed by preprocessing.

A simple and high accuracy model for time-history evolution of excessive pressure distribution of magma is not yet known. Thus, assuming that occurrence of abrupt change is rare, we use the trend model

$$\begin{aligned} \frac{d\varvec{X}}{dt}= \varvec{X}+\varvec{v}_t, \end{aligned}$$

(10)

in this paper. Here, $\varvec{X}(t)=(c_1(t),\cdots ,c_M(t))$ and $\varvec{v}_t$ is system error following the Gauss distribution. Temporal discretization of Eq. (10) with the finite-difference method leads to

(11)

Here, it indicates the time step number. The observation equation becomes

(12)

Using the scaling parameter $\alpha ^2$, the covariance of the evolution equation can be expressed as

$$\begin{aligned} \varvec{Q}_{it}=\alpha ^2\varDelta t\varvec{I}, \end{aligned}$$

(13)

where $\varvec{I}$ is the identity matrix and $\varDelta t=t_{it}-t_{it-1}$. The covariance of observation noise becomes

$$\begin{aligned} \varvec{R}_{it}=\left( \begin{array}{cc} \sigma _1^2\varvec{\varSigma }_1&{}\varvec{0}\\ \varvec{0}&{}\sigma _2^2\varvec{\varSigma }_2 \end{array}\right) , \end{aligned}$$

(14)

where $\sigma _1^2,\varvec{\varSigma }_1$ and $\sigma _2^2,\varvec{\varSigma }_2$ are the covariance and correlation matrices of GNSS and InSAR observation errors, respectively. As the state-space Eqs. (10) and (11) are expressed as linear models following the Gauss distribution, we can use the Kalman filter algorithm for data assimilation. We use a Kalman smoother with length 200.

In order to conduct data assimilation of the above, computations of Green’s functions for surface ground displacement response for each basis function is required. We use the finite-element method for this computation in this study.

2.2 Computation of Green’s Functions Using the Finite-Element Method

We target the linear elastic crustal deformation under volcano expansion. The governing equation is

$$\begin{aligned} \sigma _{ij,j}+f_i=0, \end{aligned}$$

(15)

where

$$\begin{aligned} \sigma _{ij}=\lambda \epsilon _{kk}\delta _{ij}+2\mu \epsilon _{ij}. \end{aligned}$$

(16)

Here, $\sigma ,f$ are stress and external force, $\lambda$, $\epsilon$, $\delta$, $\mu$ are the first Lame coefficient, strain, the Kronecker delta and the shear modulus, respectively. By spatial discretization with the finite-element method, the governing equation becomes a linear system of equations

$$\begin{aligned} \varvec{K}\varvec{u}=\varvec{f}, \end{aligned}$$

(17)

where $\varvec{K},\varvec{u},\varvec{f}$ are the global stiffness matrix, nodal displacement, and external force vectors, respectively.

We use the method in [2] for expressing external force, where an equivalent displacement field for excessive pressure in a void in the magma chamber is expressed using infinitesimal spheroidal pressure sources. This method expresses body force as a sum of three perpendicular dipole moments. For an expansion source with dipole moment of $f_0$ at coordinate $\varvec{x}'$, the body force $\varvec{f}_c$ can be expressed as

$$\begin{aligned} \varvec{f}_c=f_0\nabla _{\varvec{x}}\delta (\varvec{x}-\varvec{x}'), \end{aligned}$$

(18)

where $\delta (\varvec{x})$ is the Dirac’s delta function. The sizes of the dipole moments are given by Bonafede and Ferrari [1] as

$$\begin{aligned} f_0=a^3\varDelta P\frac{\lambda (\varvec{x}')+\mu (\varvec{x}')}{\mu (\varvec{x}')}\pi , \end{aligned}$$

(19)

where a and $\varDelta P$ are the radius and the pressure increment of the spherical pressure source. The external force vector in the finite-element analysis becomes

$$\begin{aligned} \varvec{f}=\sum _e\int _{\varOmega _e}[\varvec{N}]_e^T\varvec{f}_cd\varOmega , \end{aligned}$$

(20)

where $\varvec{N}$ and $\varOmega$ are the shape functions and element domains. By applying Eq. (18) to Eq. (20), the external force vector becomes

$$\begin{aligned} \varvec{f}=\sum _ef_0\nabla [\varvec{N}(\varvec{x}')]_e. \end{aligned}$$

(21)

The outer force for a basis function for spatial distribution $\varvec{B}_k(\varvec{x})$ becomes

$$\begin{aligned} \varvec{f}_k=\int \sum _e\varvec{B}_k(\varvec{x})\nabla [\varvec{N}_k(\varvec{x})]_edV. \end{aligned}$$

(22)

By using this as the right-hand side vector of Eq. (17) and solving the linear system of equations, we can obtain $\mathbf{u}$ which corresponds to the Green’s functions.

Most of the computation cost involved in finite-element analysis is generation of the finite-element model and the solver cost. We use the method of [8] to robustly generate finite-element mesh from digital elevation maps of crust structures. In our method, the spherical pressure sources are expressed without explicit modeling of the void regions in magma chambers. Thus, we do not require mesh that reflects the magma chamber geometry; we can use a single mesh reflecting the crustal structure for analyzing any pressure source. On the other hand, we need to solve the linear equation for each basis function. Thus, we can expect that the number of solver computations becomes more than the number of models to be generated. Thus, reduction of the time required for solving each of the linear set of equations becomes important for reduction of the cost of the whole procedure.

We use the method developed by Yamaguchi et al. [14] for the finite-element solver. This solver was originally developed for island-scale crustal deformation analysis for a given inter-plate fault-slip on CPU based large-scale computing environments, and was ported to GPU environments using OpenACC. The solver is based on the Adaptive Conjugate Gradient method, where the preconditioning equation is solved roughly using another Conjugate Gradient solver [6]. Here we call the iterations of the original solver as the outer loop, and the iterations of the solver used in solving the preconditioning equation as the inner loop. As the inner loop is only used as a preconditioner, we can solve it roughly. Thus, we can implement methods that can reduce computational cost in the inner loop. In this solver, we used mixed precision arithmetic and the multi-grid method. While the outer loop is computed in double precision, the inner loop is computed using single precision. By assigning suitable thresholds for the inner loop solvers, we can move most of the computation to the low-precision arithmetic parts which required less computational cost. The geometric multi-grid method and the algebraic multi-grid method are used for the multi-grid preconditioning. These methods lead to reduction in computation and data transfer; thus, it is effective for reducing computational cost on both CPUs and GPUs.

The sparse matrix-vector product $\varvec{K}\varvec{u}$ is the most time-consuming kernel in the preconditioned conjugate gradient method. In the solver, we use the Element-by-Element (EBE) method for the sparse matrix vector products except algebraic multi-grid. The multiplication of matrix ${\varvec{K}}$ and vector $\varvec{u}$ is obtained by adding up element-wise matrix-vector products as

$$\begin{aligned} \varvec{f}\leftarrow \sum _e\varvec{Q}_e\varvec{K}_e\varvec{Q}^{T}_e\varvec{u}, \end{aligned}$$

(23)

where $\varvec{f}$ is the resulting vector, $\varvec{Q}_e$ is the mapping matrix from the local node number to the global node number, and $\varvec{K}_e$ is the e-th element stiffness matrix. By computing the element stiffness on the fly instead of storing it in memory, we can reduce the amount of memory access, and thus improve computational performance. As the matrix for algebraic multi-grid is generated algebraically, the EBE method cannot be applied. Thus, we use $3 \times 3$ block compressed row storage for this part. Although $\varvec{Q}_e$ is shown in matrix form in Eq. (23), it corresponds to random access in the actual code. Random read access to memory is involved for the right-hand side vector $\varvec{u}$ part, and random write access to memory is involved for the left-hand side vector $\varvec{f}$ part. Thus, the EBE method consists of element multiplications and random memory access. Generally, the performance on GPUs deteriorate with random memory access. Thus, in order to reduce the randomness of data access, we solve several equations at the same time. If the matrix is common, the EBE computations for several vectors can be expressed as

$$\begin{aligned}{}[\varvec{f}^1,...,\varvec{f}^n]\leftarrow \sum _e\varvec{Q}_e\varvec{K}_e\varvec{Q}^{T}_e[\varvec{u}^1,...,\varvec{u}^n]. \end{aligned}$$

(24)

Here, n vectors are solved at the same time. Randomness of memory access is reduced as the read/write are conducted in consecutive address with the length of the number of vectors. This is expected to remove the performance bottleneck of the EBE kernel.

The approach of solving several vectors at the same time was originally developed for solving the crustal deformation to inter-plate fault-slip in a short time on GPUs. Although the target problem of this paper is different, it is the same from the viewpoint that displacement response to several sources are solved. Both of the discretized equations have the same stiffness matrix as the sources are only reflected to the external force vectors. Thus, we can expect applying the multiple vector method to our problem. Thus, we extended the multiple vector method such that point sources can be computed for the volcano expansion problem. By computing multiple Green’s functions for multiple point sources, we reduce the time required for computation of each Green’s function.

3 Verification

Using the developed method, we check the convergence of numerical solution and performance of the solver. Here we use an IBM Power System AC922 with 2 16-core IBM POWER9 2.60 GHz CPUs and 4 NVIDIA Tesla V100 GPUs. We use MPI to use all GPUs during computation, and use MPI/OpenMP to use all cores in the case of CPU computation.

We compared the numerical solution with Mogi’s analytical solution [10] corresponding to response of elastic half space for a spherical pressure source. The compute domain is , , , with the input source at (0 km, 0 km, $-4$ km). The dipole moment size is $1.0\times 10^{16}\,\mathrm {Nm}$, with crust material properties of $V_p=5.0\,\mathrm {km/s}$, $V_s=2.9\,\mathrm {km/s}$, and density of $2.60\,\mathrm {g/cm^3}$. The finite-element model is generated with resolution of 500 m. From Fig. 1, we can see that the obtained solution follows the analytical solution, with relative error less than 0.02 at coordinate (0, 0, 0) km. We can see that the analysis method is converged to the correct solution.

Fig. 1.

Comparison of analytical solution (solid line) and finite-element results (circles) at horizontal and vertical directions at coordinate (x, 0, 0).

Next, we measure performance for solving 12 point sources using GPUs. For comparison, we measured elapsed time for solving 1 vector at a time on GPUs, and 1 vector at a time on CPUs. Fig. 2 shows the elapsed time per vector. We can see that with the use of GPUs and the method for computing multiple vectors at the same time, we can accelerate computation by 35.5-fold from the CPU version. The elapsed time for the developed method was 23 s per vector, which corresponds to $23\,\mathrm{s} \times 10^4= 2.66$ days for computing $10^4$ equations. We can see that many case analyses can be conducted in practical time by use of the developed method.

Fig. 2.

Performance comparison of the elapsed time for solving 1 vector

4 Data Assimilation of Magma Source of Shinmoedake

We conducted data assimilation of Shinmoedake mountain using the developed method. Shinmoedake is an active volcano in the Kirishima-mountain range in Kyushu region of Japan. Recently, supply of magma and contraction due to eruption is repeated, and thus the state of magma chamber is changing. Thus, the estimation of the state of magma source is of high demand. In this section, we generate an artificial observation data of a magma eruption process, and use it to confirm that data assimilation of magma source of Shinmoedake is possible. The target region is of size $-50$ km $\le$ x $\le$ 50 km, $-50$ km $\le$ y $\le$ 50 km, $-100$ km $\le$ z $\le$ 1.4 km. Coordinate $(x, y, z) = (0, 0, 0)$ is set to the position of Shinmoedake according to the map of Geospatial Information Authority of Japan [5] with elevation set to the reference ellipsoid of the World Geodetic System (GRS80). We use the underground structure geometry and material properties given in the Japan Integrated Velocity Structure Model by the Headquarters for Earthquake Research Promotion [7]. The target region consists of 14 layers with material properties in Fig. 3. Figure 4 shows the generated finite-element model with minimum element size of $ds=500$ m. The total degrees-of-freedom was 4,516,032 and the total number of second-order tetrahedral elements was 2,831,842.

We express the distribution of magma source $\varvec{B}_k(\varvec{x})$ using a 12 noded interpolation function shown in Fig. 5. The coordinates of nodes are given in Fig. 6. Nodes 5, 6, 7, and 8 models the volcanic vent, and the magma distribution in these sections are interpolated with one-dimensional linear functions. The other nodes at plane $z=-6$ km model the magma chamber; these sections are interpolated with two-dimensional linear functions. Each basis function is normalized such that the integration $\int \varvec{B}_k(\varvec{x}) d\varvec{x}$ over the domain becomes 1 N m. Horizontal resolution is set based on resolution of GNSS station. We assume a hypothetical trend in the magma source, and use this as the target for estimation. The target trend is shown in Fig. 8, where only nodes 2, 5, 6, 7, and 8 have non-zero values. This trend is generated such that the magma pressure at the bottom of the vent is increasing, and the magma is ascending through the vent. Considering that displacement of 4 cm is observed in EBINO-MAKIZONO GNSS baseline length in one year in 2011 [9], we set the trend such that the maximum displacement becomes the order of 10 cm. We use time step increment of $\varDelta t=0.25$ day, and total number of time steps as $T=1000$. The scaling parameter for the system noise in Eq. (13) is set to $\alpha ^2=4\times 10^{32}\,\mathrm {N}^2 \mathrm {m}^2/\mathrm {day}$. This value is set such that it becomes similar to the value of covariance at the abrupt change in the target trend.

We obtain the observation data by applying the Green’s function to the magma source. We assume that GNSS data is available every 6 h at 8 observation points located near Shinmoedake (GEONET GPS-based Control Stations of Geospatial Information Authority of Japan [4]), and assume that InSAR data is available every 14 days at $60 \times 60$ km region with 500 m mesh resolution (total of 2001 observation channels). We assume that GNSS data is available every 6 h at 8 observation points located near Shinmoedake (GEONET GPS-based Control Stations of Geospatial Information Authority of Japan [4]). The coordinates of GNSS reference stations are given in Fig. 7. We also assume that InSAR data is available every 14 days at $60 \times 60$ km region with 500 m mesh resolution (total of 2001 observation channels). The observation error is assumed to be $\sigma _1=1\,\mathrm {mm^2}$, $\sigma _2=5\,\mathrm {mm^2}$. The correlation matrix $\varvec{\varSigma }_1$, $\varvec{\varSigma }_2$ are assumed to be identity matrices with noise applied to the artificial observation data. We use the same observation data and correlation matrix in the data assimilation process. We use $\theta =30^\circ$ and $\phi =10^\circ$ assuming a typical polar orbit satellite for the InSAR observation directions. Although InSAR data is known to involve correlation between data, we neglect this correlation for simplicity. In this section, we conduct estimation for three cases with different observation data: case (a) using GNSS, case (b) using InSAR, and case (c) using both GNSS and InSAR.

Fig. 3.

Material properties of layers in the target domain. $V_p$ is P wave velocity, $V_s$ is S wave velocity, and $\rho$ is density.

Fig. 4.

Overview of the finite element model.

Fig. 5.

Location of the input points. The magma distribution is interpolated with one-dimensional or two-dimensional linear functions, and input points represent control point of the magma distribution.

Fig. 6.

Location of nodes in km shown in Fig. 5

Figure 9 shows the estimated results for each case. While node 2 representing the magma chamber is estimated for all cases, estimation at nodes 5 and 6 corresponding to the vent differs among the cases. We can see that the movement of magma in the vertical direction is less sensitive compared to its horizontal movement, and thus estimation of vertical distribution is more challenging. In order to compare the results quantitatively, we use the root mean square error $\mathrm {RMSE}=\sqrt{\frac{1}{TM}\sum _{t=1}^T\Vert {\mathbf{X}}_{t}-\hat{\mathbf{X}}_{t}\Vert ^2},$ where $\hat{{\mathbf{X}}}_{1:T}$ is the estimation obtained by data assimilation, and M is the number of basis functions ($M=12$). From values in Fig. 9, we can see that RMSE is the smallest for case (c). From this, we can see that the accuracy of inner state estimation is improved with combination of available data when compared to using a single observation source.

In order to discuss the effect of the crust structure, we finally conducted data assimilation for the case disregarding the three-dimensional crust structure. A finite-element model with flat surface and homogeneous material properties is generated, and the same data assimilation procedure was conducted. Here, material number 14 in Fig. 3 was used, and the surface elevation was set to the mean elevation of the domain. The estimation results using both GNSS and InSAR is shown in Fig. 9 (d). Even if both of the observation data is used, the estimated trend was completely different from the input trend when disregarding the crust structure. In addition, the estimation results vary with the observation data; thus, periodic peaks appear when InSAR data is included. Note that the homogeneous Green’s function works for qualitative estimation of single source, as described in [11]. However, the error in Fig. 9 (d) was significantly larger than the case considering the crustal structure; we can see that consideration of the crust structure is important for high resolution estimation of the target problem. From here, we can see that the use of Green’s functions computed by finite-element modeling proposed in this paper is important.

Fig. 7.

Coordinates of GNSS reference stations in km. z is set to the surface level.

Fig. 8.

Target input trend in strength of each dipole. Nodes that are not indicated in the figure are set to have zero dipoles.

Fig. 9.

Estimated strength of each dipole (N m) with the observation data (a) GNSS alone (b) InSAR alone (c) GNSS & InSAR on the three-dimensional heterogeneous models. (d) is the result using Half space model with GNSS & InSAR. The numbers in the legend correspond to nodes shown in the Fig. 6. A result closer to the input trend in Fig. 8; a lower RMSE is preferable.

5 Closing Remarks

In this paper, we conducted basic study towards magma state estimation using time-history crust deformation data observed at the surface. As the crust involves three-dimensional heterogeneity and that the material properties of the crust changes with the state of magma, data assimilation using many-case three-dimensional finite-element analysis is expected to be suitable. In this study, we developed a fast GPU-accelerated solver for computing crustal deformation for excessive pressure sources of the magma. The elapsed time for solving a single Green’s function was 23 s; corresponding to solving $10^4$ equations in 2.6 days. We can see that many case finite-element analysis has become possible within a reasonable time frame. In the application example, we showed that the magma source trend in a three-dimensional heterogeneous crust structure can be estimated using an artificial observation data set. Data assimilation based on a homogeneous half space model resulted in wrong results; thus, we can see the importance of considering the heterogeneous crust structure and the use of finite-element method capable of modeling three-dimensional heterogeneity. As the data assimilation results were shown to differ with the crust structure, we plan to incorporate subsurface structure and topography with much more heterogeneity than the stratified structure used in Sect. 4, and plan to estimate the effect of uncertainty in the underground crust structure in our future work. Furthermore, we plan to extend the method for application to problems with material nonlinearity of the crust.