Application of the inversion method to a GPS network for estimating the stress increment in Japan

Summary

This paper examines the applicability of an inversion method of identifying stress distribution to the Japanese Islands, whose deformation is measured by the nationwide GPS array. The present inversion method is different from ordinary inversion, which requires a parametrization of constitutive relations, even though it is applicable only in the 2-D state. It is rigorously shown that stress can be found without using constitutive relations, and the validity and usefulness are verified by carrying out a numerical simulation and a model experiment for a certain class of materials. For the GPS data, the inversion method can predict the stress increment distribution that is associated with the measured displacement increment. Under assumptions of a plane stress state for the incremental deformation and no volumetric inelastic deformation, the stress increment is computed from the GPS array data measured during 1997 and 1998; from the viewpoint of numerical analysis, the validity of the solutions is examined by changing the discretization scale, boundary conditions and other parameters. Regional constitutive relations are deduced from relations between the measured strain increment and the inverted stress increment. While further studies are definitely needed, the results of the numerical simulation suggest the potential usefulness of the present inversion method. The limitations of the inversion method when applied to the Japanese Islands are discussed, and necessary modifications are mentioned.

constitutive relations, GPS, inverse analysis, inverse problem, stress estimation

1 Introduction

The Geographical Survey Institute of Japan (GSI) has been operating a nationwide Global Positioning System (GPS) array since 1997 (see e.g. Kato et al. 1998). The primary task of this array is to monitor the regional deformation in the Japanese Islands. Continuous monitoring of crustal movement at such a high spatial resolution will provide valuable information for predicting the occurrence of huge earthquakes if anomalous behaviour is observed in some regions. For earthquake prediction, however, monitoring of the stress state is more useful since earthquakes are failure phenomena that are generated by stress. The data measured by the GPS network will be used for the simulation of crustal deformation as well. For reliable simulation, a key issue is accurate modelling of the crustal structure, and hence regional constitutive relations should be estimated by using the GPS network.

Identification of regional stress and constitutive relations is a typical inverse problem that is ill-posed. Data available for the inversion are limited compared with unknowns. It is certainly possible to formulate the inverse problem—the GPS data are used to identify the most suitable stress distribution and its increment and to find the parameters of the constitutive relations if some form of the relations is assumed. However, such inverse analysis is not easy since a huge number of unknowns will have to be identified and relations between these unknowns and the measured data will be quite complicated.

Identification of the stress state and constitutive relations is a popular research topic in the field of applied mechanics (see e.g. Tanaka & Bui 1992; Bui 1995; Tanaka & Dulikravich 1998). Recently, a new inversion method was proposed by Hori & Kameda (1998) and Hori et al. (1995). While it is applicable only to a body in a 2-D plane stress or strain state, the inversion method can find a stress distribution using a strain distribution even if relations between strain and stress are not known. The local constitutive relations, linear or non-linear, are then determined by examining relations between the measured strain and the inverted stress. Identifying stress from strain without constitutive relations sounds strange. However, the key point of the inversion method is the fact that three in-plane stress components satisfy two equations of equilibrium. Hence, if one condition is found, three components can be determined.

The inversion method was originally aimed at identifying local constitutive relations near cracks or defects that were initiated in a small material sample. This method is principally applicable to the Japanese Islands, if the deformation state is approximated as a plane stress state during a short period when the GPS array measures the crustal deformation. If the inversion method works, the stress increment associated with the measured strain increment may be identified. Furthermore, regional constitutive relations may be determined from relations between the measured strain increment and the inverted stress increment even though they are for a 2-D state. Such regional constitutive relations help us to construct a reliable model of the Japanese Islands (see Fig. 1).

Figure 1

Inversion of stress increment distribution and constitutive relations.

The objective of this paper is to examine the applicability of the new inversion method to the Japanese Islands. It is expected that some modifications to the inversion method will be needed in order that is can be applied to the GPS array data, since the inversion method was originally developed for a sample of materials whose properties were known to some extent. The content of this paper is as follows. In Sections 2 and 3, the new inversion method is formulated and its usefulness is demonstrated; the results of a numerical simulation and a model experiment are presented in Section 3. In Section 4, we apply the inversion method to the GPS array data. The stress increment is computed for the strain increment that was measured during 1997 and 1998. It is shown that the results obtained support the basic applicability of the new inversion method to the Japanese Islands, although several limitations are found. Some concluding remarks are made in Section 5.

2 Formulation Of The Inversion Method

First, we briefly summarize the new inversion method (see Hori & Kameda 1998 and Hori et al. 1995). We consider a thin body, V, which is in a plane stress state. The surface and its boundary are denoted by S and ∂S, respectively, and the x₁- and x₂-coordinates are taken on S. The problem setting is as follows. The body consists of an unknown non-linear elasto-plastic material, and displacement and traction, denoted by u_i and t_i (i=1, 2), are measured on S and along ∂S, respectively (see Fig. 2).

Figure 2

Body in a plane stress state.

The constitutive relations of the elasto-plastic material are expressed by using the flow rule, the decomposition of the strain increment, ε˙_ij, into an elastic part and an inelastic part, 

1

where superscripts e and p stand for the elastic and plastic parts. The elastic and plastic constitutive relations are as follows: (1) the elastic strain increment is related to the stress increment through the elasticity tensor; and (2) the plastic strain rate is given as the gradient of a yield function, f, which is a function of stress, i.e.

where λ is the plastic coefficient. While each material has its own yield function, the plastic strain increment becomes incompressible for metals if the plasticity is induced by sliding of the crystal lattice. For the present problem, therefore, we can assume that ε˙^p₁₁+ε˙^p₂₂=0.

Now, we consider three in-plane stress components, σ₁₁, σ₂₂ and σ₁₂. Since they satisfy two equations of equilibrium, their increments also satisfy 

2

 

3

Taking advantage of eq. (1), we rewrite eqs (2) and (3) using the plastic strain increment. For instance, the first term of eq. (2) becomes

where, for simplicity, c_11klε_kl and −c_11klε_kl are replaced by σ^a₁₁ and σ*₁₁. They are called the apparent stress and the eigenstress. [Eigenstress or polarized stress is often used in applied mechanics. The origin of this term is the equivalent inclusion method proposed by Eshelby (1957), who analytically solved an inclusion problem of an ellipsoidal inclusion embedded in an unbounded body subjected to far-field loading (see also Nemat-Nasser & Hori 1998).] Note that the apparent stress increment is automatically determined if the strain increment is given, i.e. 

4

Assuming elastic isotropy, we can rewrite the assumption that ε˙^p₁₁+ε˙^p₂₂=0 as 

5

Hence, eqs (2) and (3) yield 

6

 

7

Since stress is the sum of the apparent stress and the eigenstress, traction along the boundary is given as 

8

 

9

where n_i is the unit normal. Eqs (6) and (7) together with eqs (8) and (9) form a boundary value problem for σ˙*₁₁=−σ˙*₂₂ and σ˙*₁₂.

The above boundary value problem is easily decoupled. For instance, by differentiating eqs (6) and (7) with respect to x₁ and x₂, respectively, and subtracting the resulting equations, we obtain a governing equation for σ˙*₁₁. Also, we can derive boundary conditions for σ˙*₁₁, multiplying eqs (8) and (9) with n₁ and n₂, subtracting the resulting equations and using n₁²+n₂²=1. The boundary value problems for σ˙*₁₁=−σ˙*₂₂ and σ˙*₁₂ are 

10

and 

11

These linear boundary value problems have unique solutions that can easily be obtained by using numerical methods such as a finite difference method or a finite element method. Using the solutions, we can compute the increment of the three stress components as follows: 

12

3 Verification Of The Inversion Method

In order to verify the validity of the new inversion method, we carry out a numerical simulation and a model experiment. The results of the inversion are briefly summarized in this section; see Hori & Kameda 1998 and Hori et al. 1995 for a more detailed explanation.

3.1 Numerical simulation

A stress field in a rectangular material sample is studied with a numerical simulation. The elasticity is isotropic, with a Young's modulus E and a Poisson's ratio ν and the plasticity is given by a von-Mises yield function, f=τ−τ_y, with τ and τ_y being the maximum shear stress and the yield stress; see Table 1 for the values of E, ν and τ_y. The boundary conditions are as follows: (1) a vertical load is applied to part of the upper edge; (2) the bottom edge is fixed; and (3) the right and left edges are traction-free (see Fig. 3). Displacement and stress are first computed by using a finite element method, and the computed displacement and traction are used as input data for the inversion. Nodes of the finite element method analysis are regarded as the points where these data are given. It should be emphasized that no information is provided about the plastic properties, although the elastic properties (E and ν) are known.

Table 1

Properties of the material sample.

Figure 3

Numerical simulation of a material sample.

For each increment of the applied load, eqs (10) and (11) are solved, and the stress increment is computed from eq. (12). A typical example of the stress increment inversion is shown in Fig. 4, where (a) is the exact distribution of σ˙₁₂ and (b) is the distribution obtained from the inversion. Both graphs show the stress concentration near the points where the vertical load is applied, and the agreement of these distributions is satisfactory. Since stress and strain are given, we seek the plastic constitutive relations and plot the plastic strain increment in the principal stress plane. (In this manner, we can visualize the yield function f as a surface on which the plastic strain rate is located.) The plot is shown in Fig. 5, in which (a) is the exact plastic strain increment, showing an ellipse that corresponds to the von-Mises yield surface, and (b) is the computed inelastic strain rate. Although the ellipse in Fig. 5(b) is less clear, we can evaluate the yield stress using an edge of the ellipse and obtain τ_y=0.001±0.0002 MPa, which is close to the exact value of 0.001 MPa.

Figure 4

Comparison of stress distribution. (a) Exact stress increment (dσ₁₂). (b) Predicted stress increment (dσ₁₂).

Figure 5

Comparison of yield locus. (a) Exact distribution of plastic strain increment. (b) Predicted distribution of plastic strain increment.

3.2 Model experiment

Next, we carry out a model experiment with a gelatin sample (see Fig. 6). The sample is disc-shaped, with a radius of 30 cm and a thickness of 1 cm. The torsional loading is applied at the bottom by rotating a coaxial circular part of radius 15 cm, with the peripheral site being fixed. Displacement on the surface is measured by using image analysis. Markers 0.5 cm apart forming a grid pattern are put on a rectangular region 7 cm by 8 cm on the surface, and this region is recorded by a video camera. For each image frame of the video records, the location of each marker is identified and displacement is computed. As a typical example, a video image and a maximum shear strain distribution are plotted in Figs 7 and 8. The formation of shear cracks and the strain concentration are observed. Traction along the boundary of the rectangular domain is set to zero, since the overall applied load does not increase significantly during the formation of cracks. The elastic properties are measured from a simple shear test.

Figure 6

Schematic view of model experiment.

Figure 7

Example of video images.

Figure 8

Maximum shear strain distribution obtained by image analysis.

Our concern is the local constitutive relations, and the stress–strain relations near the crack are plotted in Fig. 9, in which (a) and (b) show the normal and shear components (σ₁₁–ε₁₁ and σ₁₂–ε₁₂ relations), respectively, where the x₁- and x₂-axes are taken as the horizontal and vertical directions. While it appears that the relations between these components differ from place to place, the stress and strain invariants show some common relations. The maximum shear stress is plotted with respect to the maximum shear strain in Fig. 10, in which A–E are the five elements shown in Fig. 7. For elements A and B, through which the crack passes, the maximum shear strain increases even though the maximum stress attains an almost constant value. For the other elements, the maximum shear stress increases almost linearly. This is because the principal stress is compressional and no failure takes place. The fact that the constitutive relations common to these five elements are obtained supports the validity of the proposed inversion method.

Figure 9

Stress–strain relations in elements. (a) Normal stress–strain relation. (b) Shear stress–strain relation.

Figure 10

Relation between maximum shear stress and strain.

4 Application Of The Inversion Method To The Japanese Islands

The results of the numerical simulation and the model experiment verify the validity and the usefulness of the new inversion method formulated in Section 2. We now examine the applicability of the inversion method to the Japanese Islands, using the displacement increment (displacement rate or velocity), which is measured by the nationwide GPS array. It should be emphasized that the new inversion method is applicable only to a 2-D plane stress state. Although the stress state of the Japanese Islands varies vertically due to gravity effects and crustal structure, we assume a plane stress state for the deformation increment, at least for the shallow part of the islands, during a short period measured by the GPS array. Therefore, the inversion method provides the distribution of the stress increment that corresponds to the horizontal displacement increment measured by the GPS.

In the present analysis, the displacement data measured by the GPS array were modified. The least-squares prediction was applied to get rid of errors inherent to each GPS station. This prediction is to find the most suitable correlation length for the displacement increment of all GPS stations (see El-Fiky & Kato 1999a,b). The distribution of the displacement increment is interpolated by taking the sum of the contributions from all stations. In this manner, we can obtain a spatially smooth distribution of displacement and strain, although some irregular regional change in displacement or strain increment may be underestimated. The original GPS data are given in GSI (2000)

Attention must be paid to the assumption made in the present inversion method, i.e. eq. (5). While this assumption comes from the incompressibility of the plastic strain, a similar incompressibility can hold for a body that has sliding cracks. To show this, we consider a 2-D domain, D, with a boundary ∂D, and we denote by &OHgr; a set of cracks located in D. When the cracks have a displacement gap [u_i], according to the Gauss theorem the average strain over D is given as

where sym stands for the symmetric part. The first term on the right-hand side is the apparent strain of the domain, and the second term is the inelastic strain that is induced by crack sliding (see Nemat-Nasser & Hori 1998). When D is uniformly elastic and the elasticity tensor is c_ijkl, the average stress of this domain, denoted by &Sgr;_ij, is given as c_ijklε¯_kl, hence we have

where ε_ij and ε^c_ij are the apparent strain and the inelastic strain, respectively. This relation clearly shows the more ductile behaviour of D due to the presence of cracks. Furthermore, when the cracks slide, the displacement gap satisfies n_i[u_i]=0, and hence the inelastic strain must satisfy

If the elasticity c_ijkl is isotropic, this incompressibility of ε^c_ij yields eq. (5) when σ*_ij is given as σ*_ij=−c_ijklε^c_kl. Therefore, we can assume that eq. (5) holds, provided that (1) the regional inelastic deformation is mainly due to the fault movement and (2) the Japanese Islands have more or less uniform and isotropic elastic properties. It should be mentioned that in the current example a Young's modulus and Poisson's ratio of E=1 MPa and ν=0.3 are used. This value of E is used as a reference, hence the absolute value of local stress is not meaningful.

4.1 Preliminary check of applicability

In order to examine the applicability of the new inversion method to the Japanese Islands, we make three preliminary checks. First, we examine the convergence of the numerical solutions. A finite element method is employed to solve eqs (10) and (11); square isoparametric elements with nine nodes are used. Since a continuous distribution of strain is given, we can change the element size freely. In Figs 11 and 12, the distributions of the first and second stress invariants, 

13

Figure 11

Convergence of the distribution of σ˙.

Figure 12

Convergence of the distribution of τ˙.

are plotted for three cases of discretization. Figs 11 and 12(a), (b) and (c) show the results for elements of side length &Dgr;=0.5, 0.25 and 0.125°. The difference between the discretizations for 0.25 and 0.125° is barely discernible. Similar results are obtained for three stress components (σ₁₁, σ₂₂ and σ₁₂). Therefore, the convergence of the numerical solutions is verified.

Next, we examine the effects of the boundary conditions. Due to the linearity, the boundary value problem for, say, σ˙*₁₁, is decomposed into two parts, i.e. σ˙*₁₁=σ˙*^(a)₁₁+σ˙*^(t)₁₁, where σ˙*^(a)₁₁ and σ˙*^(t)₁₁ satisfy

and

σ˙ *^(a)₁₁ corresponds to the apparent stress on the surface and σ˙*^(t)₁₁ to the traction along the boundary. They give the stress increment as σ˙_ij=(σ˙*^(a)_ij+σ˙^a_ij)+σ˙*^(t)_ij, with σ˙*^(a,t)₂₂=−σ˙*^(a,t)₁₁. Since the traction increment along the Japanese Islands is not known, we consider σ˙*^(a)_ij only. Choosing two configurations for S, we solve the boundary value problem for σ˙*^(a)₁₁=−σ˙*^(a)₂₂ and σ˙*^(a)₁₂. Figs 13, 14 and 15 plot the three components of the stress increment, σ˙₁₁, σ˙₂₂ and σ˙₁₂, respectively; (a) and (b) are for larger and smaller configurations, denoted as Type 1 and Type 2, respectively, which are indicated by a solid line in the figures. Except for regions close to the boundary, the distributions are similar to each other. It should be noted that if the traction increment is uniform, that is, ⃛_i is given as ⃛_i=n_jσ˙′_ji with constant σ˙′_ji, then the corresponding eigenstress σ˙*^(t)_ij becomes uniform. Under the assumptions that the actual traction increment along the Japanese Islands is more or less uniform and that the apparent stress increment corresponds to the regional strain increment, the non-uniform distribution of the stress increment will be mainly due to the regional strain increment. This may explain the lower dependence of the inverted stress distribution on the boundary conditions, as shown in Figs 13–15.

Figure 13

Effects of boundary conditions on the distribution of σ˙₁₁.

Figure 14

Effects of boundary conditions on the distribution of σ˙₂₂.

Figure 15

Effects of boundary conditions on the distribution of σ˙₁₂.

Finally, we examine the effect of the reference elasticity that is used for computing the apparent stress increment. It should be noted that the apparent stress increment changes depending on the reference elasticity. The eigenstress increment, however, compensates such a change and can produce the exact stress increment for any choice of reference elasticity, provided that the assumption eq. (5) is correct. Since we assume isotropy, Poisson's ratio ν is changed; the apparent stress is linear to Young's modulus E, and the change in E leads to a change in the scale of the stress increment, not in the distribution. We examine four cases of ν, and plot the distributions of the first and second invariants in Figs 16 and 17, in which (a), (b), (c) and (d) are for ν=0.1, 0.2, 0.3 and 0.4, respectively. While global patterns of the distribution are similar, some regional stress increments are different for the case of ν=0.1 and 0.4. For the other two cases, the distributions are almost the same but not identical. Therefore, we can see some effects from the reference elasticity. This will be discussed below.

Figure 16

Effects of reference elasticity on the distribution of σ˙.

Figure 17

Effects of reference elasticity on the distribution of τ˙.

4.2 Results of the computation and discussion

We now apply the new inversion method to compute the distribution of the stress increment. The following conditions are set for the computation: the discretization is &Dgr;=0.125, the configuration is Type 1 as shown in Figs 13 and 14, and Young's modulus and Poisson's ratio are E=1 MPa and ν=0.25. The boundary traction increment is set to zero. In Figs 18 and 19, the first and second invariants of the strain and stress increments are plotted; the strain invariants are defined as

Figure 18

Comparison of the distribution of ε˙ and σ˙.

Figure 19

Comparison of the distribution of γ˙ and τ˙.

The overall patterns of the stress increment invariants, σ˙ and τ˙, are similar to the corresponding strain increments, ε˙ and γ˙, respectively; for instance, in Fig. 18, regions where ε˙ is tensile or compressive have tensile or compressive σ˙, respectively, and in Fig. 19, regions with large γ˙ have large τ˙ as well. However, the regional distribution is slightly different from place to place, and the degrees of concentration for strain and stress increments are not the same. This clearly shows that the apparent stress increment, which is computed from the strain increment through the assumed elasticity, does not satisfy the equilibrium, and hence the eigenstress increment is generated such that the stress increment satisfies the equilibrium. Fig. 20 shows the orientation of the maximum shear strain and stress increments, φ and θ, i.e.

Figure 20

Comparison of the distribution of φ and θ.

Like in Figs 18 and 19, the regional distributions of these orientations are different, even though the overall patterns are similar to each other.

It should be mentioned that the stress increment invariants show some peculiar distributions near the boundary. This is partly because the boundary traction increment is neglected in solving the problems. Also, the strain data used in the inversion are not suitable for these regions since the strain increment is extrapolated there (see Figs 18a or c and 19a or c).

After the stress increment has been computed, we seek regional constitutive relations using relations between the measured stain increment and the inverted stress increment. It should be emphasized that these relations are automatically obtained once the inversion is made. At the present stage, the validity of the predicted relations cannot be verified at all.

In Fig. 21, the distributions of τ˙/γ˙ and θ−φ are plotted. By definition, τ˙/γ˙ gives the regional shear stiffness, and θ−φ is a measure of the regional anisotropy, since the orientations of the principal shear stress and strain axes coincide for an isotropic material. In Fig. 21(a), there is a high shear stiffness for the east coast of Honshu Island and the south coast of Shikoku Island. These two regions have small θ−φ, as shown in Fig. 21(b). These two observations mean that the east coast of Honshu Island and the south coast of Shikoku Island are more or less isotropic and have large stiffnesses. As mentioned, the validity of this result has not been verified; Fig. 21 is a demonstration of the potential usefulness of the new inversion method.

Figure 21

Distribution of regional stiffness and anisotropy.

We now consider the limitations of the present inversion method. It appears in Fig. 21(c) that some regions have quite high shear stiffnesses where τ˙/γ˙ is larger by an order of magnitude than in other regions (e.g. the east coast of Honshu Island). Since the strain increment was smoothed by the least-squares prediction (El-Fiky & Kato 1999a,b), these regions have a high strain increment. Together with the approximation of zero traction increment along the nearby boundary, this leads to an unrealistically high stress increment concentration being obtained. These two points indicate the limitations of applying the inversion method to the Japanese Islands, that is, the usage of smoothed strain and the approximation of zero traction increment. However, using the GPS data without filtering measurement noises may give rise to another difficulty in solving the boundary value problems, eqs (10) and (11). This is because the second derivatives of the apparent stress are used in the governing equations. Thus, some care must be taken in using data measured at stations of the GPS array. We are now reformulating the inversion method in a weak form such that the displacement is used as the input data for the boundary value problem (Hori & Kameda 2000), as well as developing a method to estimate the boundary conditions by considering the movement of the plate surrounding the Japanese Islands (Hori 2000).

Finally, we discuss the validity of the assumption eq. (5), which plays a key role in deriving the boundary value problems eqs (10) and (11). This assumption is strong in the sense that the bulk stiffness, σ˙/ε˙, becomes constant when c_ijkl is isotropic, i.e.

As shown in Figs 21(a) and (c), the shear stiffness given by τ˙/γ˙ varies from place to place. Also, the regional anisotropy is observed in Figs 21(b) and (d). Therefore, we have to change the condition of constant σ˙/ε˙, such that both the regional bulk and the shear stiffness vary in a harmonic manner. We can easily change the condition by regionally changing the reference elasticity. The determination of such regional reference elasticity will be treated as an optimization problem or another inverse problem of finding a suitable reference elasticity tensor. Some modification will be required in the inversion method, which was originally applied to material samples used in experiments.

5 Concluding Remarks

The results of preliminary checks and trial computation shown in the preceding section support the applicability of the new inversion method to the Japanese Islands. Since the method takes advantage of the fact that stress is in equilibrium, the distribution of the inverted stress increment is slightly different from the distribution of the measured strain increment. Although further investigation is definitely needed, the present results indicate the potential usefulness of applying the new inversion method to the Japanese Islands.

The results also show some drawbacks and limitations of the inversion method, such as the use of smoothed strain data, the treatment of the boundary conditions, and the strong assumption provided by eq. (5). To overcome these problems, we are now reformulating the inverse problem and the corresponding inversion method such that the resulting boundary value problems will be of more suitable form for application to the stress inversion of the Japanese Islands (see Hori (2000) and Hori and Kameda (2000). Comparisons with known crustal structures will be made once such modifications are completed for the inversion method.

Acknowledgments

This research was supported partially by a Grant-in-Aid for Scientific Research, the Ministry of Education, Science, Sports and Culture, Japan, and partially by the Japan Science and Technology Corporation.