Numerical techniques based on the finite difference scheme with discontinuity have been developed for problems associated with aseismic ground deformation in seismically active regions. We stress upon the applicability of such numerical techniques in solving problems in geodynamics. A long strike-slip fault is considered in a viscoelastic half space representing the lithosphere-asthenosphere system. The fault undergoes a sudden movement under the action of tectonic forces induced by mantle convection. The resulting boundary value problems have been solved with the help of numerical technique, based on a finite difference scheme with appropriate boundary conditions, developed for the purpose. The numerical techniques developed here can be modified for more general deformation problems where analytical methods become very complicated.
Small ground deformations have been observed in seismically active regions during the aseismic period in between two major seismic events. A long, vertical, surface breaking, strike-slip fault is taken to be situated in a viscoelastic half space of Maxwell type material. Mantle convection and other related phenomena, result in tectonic forces, which act on the system producing a sudden movement across the fault leading to an earthquake. Such types of models were considered in the article 1, 2, 3. They have obtained analytical solutions of displacement, stress and strain in the medium using integral transform and Green's functions techniques. In the present paper a numerical technique has been developed to study the nature of displacements, stresses and strains in the system near the fault after fault movement. The numerical technique is basically a two dimensional finite difference scheme with discontinuity along the fault plane. The advantages of numerical techniques lie in the fact that the movement across the fault as well as the nature of tectonic forces due to mantle convection can be represented more realistically following the actual observations on ground deformations.
We consider a theoretical model of the lithosphere asthenosphere system with long, vertical, surface breaking strike slip fault F of finite depth D lies in a linearly viscoelastic medium with its material of the Maxwell type.
We introduce a rectangular Cartesian coordinate system () with the free surface of the viscoelastic half space given by and axis is pointing vertically downward into the half space. axis lies on the free surface of the viscoelastic half space and extended along the upper edge (strike) of the fault F. With this choice of axes the half space is given by and the fault F by , where D be the vertical depth of the fault from free surface. A simple diagram of our model with respect to the Cartesian coordinate system is shown in the following Figure 1.
Let be the components of the displacement vector along the coordinate axes respectively and , , , , , be the components of stress while , , , , , are the components of strain. Due to long fault along axis all the stress, strain and displacement components are taken to be independent of and therefore they are functions of and t only. Under this assumption all stress, strain and displacement components separate out into two independent groups-one group consisting of ; , ; , is associated with possible strike slip movement while the other group consisting of ; ,, ; , , is associated with possible dip slip movement of the fault as explained in 4 and 5. Here we consider only the first group for studying possible strike slip movement in the medium.
We set our initial time t=0 when the entire medium is in a quasi static, aseismic state but a slow aseismic deformation continuously going on within it. We assume that two equal and opposite constant shear stresses , are acting in the medium far away from the fault plane and parallel to fault plane. These are maintained by some tectonic forces arising mainly due to mantle convection and/or other geological changes. Therefore for t 0 the displacement, stress and strain components satisfy the following relations and boundary conditions.
2.1. Stress-Strain Relations (Constitutive Equations)For the viscoelastic medium of Maxwell type the relation of stress and strain components associated only with strike -slip movement are obtained by the following relations:
(1) |
(, .
Where is the effective rigidity and ɳ is the effective viscosity of the viscoelastic medium and supposed to be constant throughout the medium.
2.2. Stress Equation of MotionAs per our consideration the entire medium is in a quasi static, aseismic state but a slow, aseismic deformation continuously going on within it. Due to this slow aseismic quasi static deformation, order of magnitude of variation of inertial forces are very small and can be neglected.
Therefore, the relevant stress equation of motion will be of the form:
(2) |
(,
From equations (1) and (2) we get,
This is satisfied if
(3) |
(, .
2.3. Boundary ConditionsOn the free surface the shear stress component must vanish and we have
(4) |
At a large vertical depth from free surface stress component remains unaltered with respect to the stress level in the medium at t=0. Therefore,
(5) |
As our earlier consideration two equal and opposite constant shear stresses , are acting in the medium far away from the fault plane and parallel to the fault plane. So we get
(6) |
.
Where, is the constant shear stress maintained by mantle convection and some other related tectonic phenomena.
2.4. Initial ConditionsLet , are the initial values (i.e. at t=0) of respectively. They are functions of and satisfy all the relations from (1) to (6) where initial time t= 0 is measured from a suitable instant when there is no seismic activity in the medium.
When there is no seismic events in the medium the displacement, stress and strain components associate with the strike slip movement are obtained by applying Laplace transform with respect to t of the constitutive equation (1), reduced equation (3), boundary conditions from (4) to (6) and finally using inverse Laplace transform as follows: 2
(7a) |
At left hand side region of the fault plane displacement, stress and strain components are given by the above relation (7a) replacing by and by
(7b) |
From the above result it is clear that the stress component is a increasing function (on the left hand side it is increasing negatively) of t. From the initial value it will ultimately reach to or as t . From the rheological behavior of the viscoelastic material near the fault F we may assume that it is capable to bear only some threshold amount of stress (say, called the critical value of stress, with ||< ||. Therefore after a certain time, say T, the accumulated stress in the medium near the fault F exceeds the critical level of stress and a fault movement will take place across the fault plane. In this model we consider a sudden movement along the fault plane due to overstep of the accumulate stress near the fault plane after the critical time T.
From the previous section we have seen that after a certain time T, the stress component , which is the main driving force in our model generated mainly due to mantle convection and some other tectonic phenomena such as gravity, internal pressure, resisting forces along the plate boundaries etc, oversteps the critical value of the stress and a dislocation (assumed to be a sudden movement) occurs in the medium across the fault F.
Therefore, all the relations from (1) to (6) are also satisfied in this case with an extra condition due to dislocation along F as follows:
(8) |
where [] is the discontinuity in displacement across the fault F, defined by [] =- (0).
Where U is the amount of dislocation on the free surface, H() is the Heaviside step function and is the depth dependence dislocation function, which should follows the following conditions for finite values of displacements, stresses and strains at any point in the medium including the point on lower edge of the fault. 2
(i) and are continuous functions of for 0.
(ii) and .
(iii) Either is continuous in 0 except for a finite number of points of finite discontinuity in 0 or, is continuous in 0 and there exist real constants m, n < 1 such that or to a finite limit as and or to a finite limit as .
Here we consider depth dependence dislocation function which satisfies all the above sufficient conditions as:
(9) |
After the critical time T the medium passes through a seismic state for a small time period, say which is very small, usually of the order of a few seconds.The seismic disturbances gradually die out and aseismic state reestablished. We study our model at the time (= T + ) immediately after restoration of new aseismic state in the medium.
The relevant boundary value problem after the fault movement can be described as follows:
The result of seismic tomography, numerical simulation of mantle convection and examination of Earth's gravitational fields suggest that horizontal movement of plates in two opposite directions can be occurred along the plate boundary due to mantle convection. On the basis of this observation, we consider two equal and opposite constant forces =constant = acting along + direction as and =constant = acting along direction as and remain unchanged after dislocation.
We are introducing finite difference method to find the displacements, stresses and strains over the region. Because of long fault along axis all the components are independent of and therefore we consider a rectangular region R ( , 0) perpendicular to the axis through origin, where and b are chosen as per our requirement. The surface breaking, long, vertical strike slip fault lies at and its dislocation extended through axis upto depth D from free surface.
In order to find the numerical solution of (3) with appropriate boundary conditions, we superimpose on R a rectangular network with equal step length h in and directions (i.e. square network). The nodal points are given by
Denote the numerical value of and at the nodal point (, by and and respectively.
For applying finite difference scheme we write equation (1) and (3) in discrete form taking equal step length h for and directions (i.e. square network), step length for time variable t and obtain the following expressions:
For equation (3) at time t,
(10) |
For first relation of equation (1),
(11) |
where and denote value of and respectively at the time t.
But, following our consideration far away from the fault (i.e. as |) stress remains constant and therefore at the boundary regions (, ) we must have =. By this condition the above expression (11) becomes
(12) |
For equation (5) as (i.e. at a large depth from free surface) at any time t,
(13) |
The horizontal displacement on the free surface (, ) at time t = only due to dislocation along the fault F is given by the following expression. 6 (by using correspondence principle)
(14) |
[].
4.1. Application of Finite Difference Method with DislocationThus, we obtain a rectangular region R(, 0) having a discontinuity in along the fault plane and satisfies the equation of motion (3) over the whole region with following conditions.
The surface breaking , vertical, long, strike slip fault F lies at the middle of the region R (i.e. , 0) and two equal opposite forces and - are acting parallel to axis at the two opposite sides and respectively of the fault. Therefore, we consider an equal and opposite (i.e. anti- symmetrical) distribution of dislocation on the two sides of the fault F.
Due to anti-symmetrical force system in the region R, we should have an anti-symmetrical distribution of dislocation on the free surface. Dislocation just at the positive side of the fault is and just at the negative side and which is consistent with formula (14). Because, from (14) as and , ) .
On the free surface total horizontal displacement is given by adding the normal displacement ( due to effect of || obtained by (7a) or (7b) with the displacement due to dislocation along the fault F after critical time T by (14). At the right hand side of R (i.e.= ) stress component =constant = ; at the left hand side of R (i.e.= ) stress component =constant = ; at the vertically downward side of R (i.e.= b) = 0. [Figure 2]
Applying the difference scheme (10) at all the internal grid points (except at the grid points along AB and CD denoted by darken triangle in Figure 2), where AB and CD are the vertical lines of the network R just to the left and to the right of the fault F respectively, and with the help of appropriate boundary conditions stated above and using relations (12), (13) at the corresponding boundary we obtain a system of linear algebraic equations.
But there is some discontinuity of , given by relation (8), due to dislocation along the fault F. Due to presence of this discontinuity along F we can't apply the difference scheme (10) at the grid points on AB and CD. We set some other linear algebraic equation at that points using depth dependence dislocation function (9) as follows:
As two equal and opposite forces are acting parallel to axis far away from the fault at two opposite sides () of the fault F which lies at the middle position (i.e. , ) of the region R (, ). Thus we obtain an anti-symmetrical force system with respect to the fault plane over the region R. Therefore we assume an anti-symmetrical distribution of dislocation at the grid points on AB and CD. We first calculate total dislocation along the fault F using relations (8) and (9), and divide this dislocation into two equal half but in opposite directions and set a new linear algebraic equation at that grid points.
Along AB:
Or equivalently we may write
Similarly, along CD:
Where is the Cartesian coordinate of the (i, j) th grid point in the region R.
By this we get a system of linear algebraic equations with the same number of equations as the number of unknowns. This system of algebraic equations can be expressed as
(15) |
For the linear partial differential equation + +++ in a regain ℜ the difference scheme will be of the form
(16) |
where A,B, ... , G are continuous functions of x and y and following 7 we get coefficient values as:
By applying the difference scheme (16) at all the nodal points in the region ℜ, we get a system of algebraic equations with the number of equations equals to the number of unknowns and this system of equations can be written as
(17) |
Under the following condition
(18) |
We have,
(19) |
i.e. off-diagonal elements of M are positive, diagonal elements are negative and M has diagonally dominance. Thus the matrix M is irreducible and by these properties we assure the convergence of system of equations (17).
In our problem the linear partial differential equation is (Laplace Equation) i.e. A = B = 1, C = E = F = G = 0 and for square network h = k. So, the condition (18) is satisfied obviously by our partial differential equation and the coefficient matrix P in equation (15) is irreducible 14 which assures the convergence of the solution of the system of equations (15).
It is important to note for the elliptic equation method of support operators (which we use for the finite difference scheme) that they automatically give us a stable finite-difference scheme. 9
For finding results using numerical computation we consider different values for our model parameters as follows:
Following the paper 2 and from the recent studies on Lithosphere-asthenosphere system we consider the rigidity (μ) and viscosity (ɳ) as
U = Dislocation across the fault on free surface = 80 cm.
For most of the major faults, faults width are found to lie in between 10 to 15 km and therefore we take
Under this critical stress value it is found that oversteps the critical value of stress after a time of 62.8 years approximately and we take T = Critical time level =62.8 years.
From equation (14) we have seen that at a horizontal distance of 150 km from the fault, in both sides of F, free surface displacement due to dislocation in becomes very small and therefore we take,
Applying finite difference scheme at each nodal point over R with previously stated modification for dislocation near the fault F, we obtain a system of linear algebraic equations . Using MATLAB (2012) we first solve this system of equations and obtain displacement at each nodal point. With the help of these displacement values we further evaluate stress at each nodal point and obtain the displacement, stress immediately after the dislocation across the fault, at a time when the aseismic state re-established in the medium. To find out the effect of dislocation on displacement and stress components we subtract displacement and stress values before dislocation (i.e. at time T) from the final results.
The Figure 3 shows variation of displacement component due to dislocation across the fault along different horizontal lines at different depth over the rectangular region R. In each cases figure is anti-symmetrical in nature with respect to dislocation line as expected. In both sides of the fault absolute value of displacement component is found to decrease as we move away from the fault and tends to zero as || > 150 km. The rate of change of displacement component with reduces in the lower region from free surface. However, near the lower edge of the fault 9km the displacement initially increases in magnitude and thereafter tends to zero as | > 150 km. At a depth of about 20 km, magnitude of the displacement component due to dislocation is found to be very small.
7.2. Displacement along Vertical LinesIn Figure 4 we show displacement component along different vertical lines parallel to the fault. In right hand side region > 0 displacement gradually decreases in magnitude along direction and ultimately tends to zero at a large depth from free surface. It is also found that the rate at which the displacement is decreasing depends upon the distance from the fault plane being highest near the fault plane. Curve corresponding to a large distance from the fault plane passes very close to axis and remains almost constant throughout the depth which signifies that far away from the fault there is a negligible effect of dislocation on displacement component and which gives a satisfactory support with similar geophysical models. Here we have plotted curves only at right hand side region > 0 but we will obtain an anti-symmetrical picture over whole region R.
In the Figure 5 we have shown variation of stress component due to dislocation along different horizontal lines at different depth over the region R and the figure is anti-symmetric as expected. In right hand side (> 0) at the upper portion of the fault, the magnitude of stress is found to be decreasing as we move away from the fault plane and from a negative value it will ultimately tends to zero as 100 km. Thus we get a stress reduction zone at the upper portion of the fault and this reduction amount decreases as we move away from the fault plane.
At the lower edge of the fault we get slightly different nature in stress. Near to the fault, stress is positive and decreases to some small negative value, and it ultimately converges to zero as we move away from the fault. Thus we obtain a stress accumulation zone near the lower end of fault. At a large depth from the free surface variation of stress due to dislocation remains unaltered as we move away from the fault.
For left hand side of the fault stress pattern can be explained in same way. But one important declaration should be noted here that in positive side of F, +ve stress means stress accumulation and -ve stress means stress reduction whereas in negative side, +ve stress means stress reduction and -ve stress means stress accumulation.
In the Figure 6 we have shown overall variation of stress due to dislocation over the region R. Where positive value means stress accumulation and negative value means stress reduction over the whole region R.
In the Figure 7 we have distinguished the whole rectangular region into two different parts; one part marked by “A” means stress accumulation region and the other part marked by “R" means stress reduction region due to dislocation along the fault plane. Therefore, for any other fault in the stress reduction region its chance of further movement reduces whereas for any other fault in stress accumulation region its chance of further movement enhances due to this dislocation across the fault F.
The objective of the present paper was to show that proper numerical technique can be developed for solving boundary value problem of complex nature having discontinuities in the variables.
The results obtained here were found to be in conformity with the similar results obtained by analytical methods.
The problems of discontinuities had been taken care of by suitable adjustment of finite difference schemes.
It is expected that such numerical methods can also be applied to more complicated types of problems related to aseismic ground deformation.
Appropriate modifications can also be incorporated in our model in case when the tectonic forces are not anti-symmetric in nature.
One of the authors (Subhash Chandra Mondal) acknowledges the library section of Department of Applied Mathematics, University of Calcutta for providing some valuable books and also to computer laboratory section for helpful cooperation with our lab-work.
The authors have no competing interests.
[1] | Debnath, P., Sen, S., “Creeping Movement across a Long Strike- Slip Fault in a Half Space of Linear Viscoelastic Material Representing the Lithosphere-Asthenosphere System”, Frontiers in Science 2014, 4(2): 21-28. | ||
In article | View Article | ||
[2] | Debnath, P., Sen, S., “A Vertical Creeping Strike Slip Fault in a Viscoelastic Half Space under the Action of Tectonic Forces Varying with Time”, IOSR Journal of Mathematics (IOSR-JM), e-ISSN:2278-5728, p-ISSN:2319-765X. Volume 11, Issue 3 Ver. I (May-Jun, 2015), pp 105-114. | ||
In article | |||
[3] | Ghosh, U., Mukhopadhyay, A., Sen, S. (1992), “On two interacting creeping vertical surface-breaking strike-slip fault in a two-layered model of the lithosphere”, Physics of the Earth and Planetary Interiors, 70, pp. 119-129 . | ||
In article | View Article | ||
[4] | Mukhopadhyay, A. et.al. (1979a), “On stress accumulating near finite fault”. Indian journal of meteorology, Hydrology and Geophysics. (Mousam), Vol. 30, pp. 347-352. (with B.B. Pal and S. Sen). | ||
In article | |||
[5] | Mukhopadhyay, A. et.al. (1979b), “On stress accumulating and fault slip in the lithosphere”. Indian journal of meteorology, Hydrology and Geophysics. (Mousam ), Vol. 30, pp. 353-358. | ||
In article | |||
[6] | Segall, P., Earthquake and Volcano Deformation, Princeton University Press, Princeton and Oxford, 1954. | ||
In article | View Article | ||
[7] | Jain, M. K., Iyengar, S. R. K., Jain, R. K. , Computational methods for Partial Differential Equation, Wiley Eastern Limited, 4835/24, Ansari Road, Daryaganj, New Delhi- 110002, 1994. | ||
In article | |||
[8] | Ghods, A.,Mir, M., “Analysis of rectangular thin plates by using finite difference method”, Journal of Novel Applied Sciences, 2014JNAS Journal-2014-3-3/260-267. | ||
In article | View Article | ||
[9] | Shashkov, M., Conservative Finite-Difference Methods on General Grids, CRC Press, Inc., 1996. | ||
In article | View Article | ||
[10] | Catlhes III, L.M., The viscoelasticity of the Earths mantle, Princeton University Press, Princeton, N.J, 1975.15 | ||
In article | |||
[11] | Chandru, M. et. al. (2017), “A Numerical Method for Solving Boundary and Interior Layers Dominated Parabolic Problem with Discontinuous Convection coefficient and Source Terms”, Differential Equations and Dynamical System, pp 1-22, | ||
In article | View Article | ||
[12] | Chift, P., Lin, J., Barcktiausen, U., “Evidence of low flexural rigidity and low viscosity lower continental crust during continental break-up in the South China Sea”, Marine and Petroleum Geology: 19, 951-970, 2002. | ||
In article | View Article | ||
[13] | Das, P., “A higher order difference method for singularly perturbed parabolic partial differential equation”, Journal of Difference Equation and Applications, Vol- 24, 2018-Issue 3, Page 452-477. | ||
In article | View Article | ||
[14] | I. Farago, Cs. Gaspar, Numerical methods to the solution of partial differential equations with hydrodynamic applications, Technical University, Budapest, 1983. [study material, not published]. | ||
In article | |||
[15] | Karato. S.I. (July, 2010), “Rheology of the Earths mantle”, A historical review of Gondwana Research, vol-18, issue-1, pp. 17-45. | ||
In article | |||
[16] | Mao, D., “A treatment of Discontinuities in Shock-Capturing Finite Difference Methods”, Journal of Computational Physics 92, 3422-445 (1991). | ||
In article | View Article | ||
[17] | Mao, D., “A treatment of Discontinuities for Finite Difference Methods”, Journal of Computational Physics 103, 359-369 (1992). | ||
In article | View Article | ||
[18] | Maruyama, T. (1966). “On two dimensional dislocation in an infinite and semi-infinite medium”. Bull. Earthquake Res. Inst., Tokyo Univ., 44, part 3, pp. 811-871. | ||
In article | |||
[19] | Rybicki, K. (1971), “The elastic residual field of a very long strike-slip fault in the presence of a discontinuity”. Bull. Seis. Soc. Am., 61, 79-92. | ||
In article | View Article | ||
Published with license by Science and Education Publishing, Copyright © 2018 Subhash Chandra Mondal, Sanjay Sen and Suma Debsarma
This work is licensed under a Creative Commons Attribution 4.0 International License. To view a copy of this license, visit https://creativecommons.org/licenses/by/4.0/
[1] | Debnath, P., Sen, S., “Creeping Movement across a Long Strike- Slip Fault in a Half Space of Linear Viscoelastic Material Representing the Lithosphere-Asthenosphere System”, Frontiers in Science 2014, 4(2): 21-28. | ||
In article | View Article | ||
[2] | Debnath, P., Sen, S., “A Vertical Creeping Strike Slip Fault in a Viscoelastic Half Space under the Action of Tectonic Forces Varying with Time”, IOSR Journal of Mathematics (IOSR-JM), e-ISSN:2278-5728, p-ISSN:2319-765X. Volume 11, Issue 3 Ver. I (May-Jun, 2015), pp 105-114. | ||
In article | |||
[3] | Ghosh, U., Mukhopadhyay, A., Sen, S. (1992), “On two interacting creeping vertical surface-breaking strike-slip fault in a two-layered model of the lithosphere”, Physics of the Earth and Planetary Interiors, 70, pp. 119-129 . | ||
In article | View Article | ||
[4] | Mukhopadhyay, A. et.al. (1979a), “On stress accumulating near finite fault”. Indian journal of meteorology, Hydrology and Geophysics. (Mousam), Vol. 30, pp. 347-352. (with B.B. Pal and S. Sen). | ||
In article | |||
[5] | Mukhopadhyay, A. et.al. (1979b), “On stress accumulating and fault slip in the lithosphere”. Indian journal of meteorology, Hydrology and Geophysics. (Mousam ), Vol. 30, pp. 353-358. | ||
In article | |||
[6] | Segall, P., Earthquake and Volcano Deformation, Princeton University Press, Princeton and Oxford, 1954. | ||
In article | View Article | ||
[7] | Jain, M. K., Iyengar, S. R. K., Jain, R. K. , Computational methods for Partial Differential Equation, Wiley Eastern Limited, 4835/24, Ansari Road, Daryaganj, New Delhi- 110002, 1994. | ||
In article | |||
[8] | Ghods, A.,Mir, M., “Analysis of rectangular thin plates by using finite difference method”, Journal of Novel Applied Sciences, 2014JNAS Journal-2014-3-3/260-267. | ||
In article | View Article | ||
[9] | Shashkov, M., Conservative Finite-Difference Methods on General Grids, CRC Press, Inc., 1996. | ||
In article | View Article | ||
[10] | Catlhes III, L.M., The viscoelasticity of the Earths mantle, Princeton University Press, Princeton, N.J, 1975.15 | ||
In article | |||
[11] | Chandru, M. et. al. (2017), “A Numerical Method for Solving Boundary and Interior Layers Dominated Parabolic Problem with Discontinuous Convection coefficient and Source Terms”, Differential Equations and Dynamical System, pp 1-22, | ||
In article | View Article | ||
[12] | Chift, P., Lin, J., Barcktiausen, U., “Evidence of low flexural rigidity and low viscosity lower continental crust during continental break-up in the South China Sea”, Marine and Petroleum Geology: 19, 951-970, 2002. | ||
In article | View Article | ||
[13] | Das, P., “A higher order difference method for singularly perturbed parabolic partial differential equation”, Journal of Difference Equation and Applications, Vol- 24, 2018-Issue 3, Page 452-477. | ||
In article | View Article | ||
[14] | I. Farago, Cs. Gaspar, Numerical methods to the solution of partial differential equations with hydrodynamic applications, Technical University, Budapest, 1983. [study material, not published]. | ||
In article | |||
[15] | Karato. S.I. (July, 2010), “Rheology of the Earths mantle”, A historical review of Gondwana Research, vol-18, issue-1, pp. 17-45. | ||
In article | |||
[16] | Mao, D., “A treatment of Discontinuities in Shock-Capturing Finite Difference Methods”, Journal of Computational Physics 92, 3422-445 (1991). | ||
In article | View Article | ||
[17] | Mao, D., “A treatment of Discontinuities for Finite Difference Methods”, Journal of Computational Physics 103, 359-369 (1992). | ||
In article | View Article | ||
[18] | Maruyama, T. (1966). “On two dimensional dislocation in an infinite and semi-infinite medium”. Bull. Earthquake Res. Inst., Tokyo Univ., 44, part 3, pp. 811-871. | ||
In article | |||
[19] | Rybicki, K. (1971), “The elastic residual field of a very long strike-slip fault in the presence of a discontinuity”. Bull. Seis. Soc. Am., 61, 79-92. | ||
In article | View Article | ||