Geometric Deformation Analysis of the Earth's Crust using GPS Observation and Non-linear Finite Element Method based on C1 Bezier Cubic Interpolation

Earth as a dynamic system is constantly undergoing the Change and variation. The changes causing factors include a wide range of dynamic processes, which consist of tidal effects, the loading on the crust and tectonic activity of tectonic plates. Accumulation of strain energy at the boundary between the plates and the surrounding area of faults are the most important consequences of movement and deformation of the tectonic plates. This energy may causes tension forces in the fracture sits, and if these forces can overcome the rocks resistance, an earthquake will occur. To measure strain energy, some mechanical tools contain strain gauges and turnkey have been installed at faults location which measure the deformations in situ. Due to the high cost of purchasing and maintenance of such these devices, based on geodetic observation, some alternative methods have been developed for study of the earth crust in local scales. Today, time measurement of the earth's crust deformations has been recognized as one of the notable branches of geodesy.
Measuring the Earth's crust deformation by geodetic methods has significant impacts on geological studies and provides a deep insight into understanding the mechanism of tectonic activities such as earthquakes and volcanoes.
Over time, measurement methods in geodesy significantly expande and space measurement methods have been replaced classic geodetic observations. The Global Positioning System (GPS) playes an important role in the development of observational methods in geodesy. By using GPS observations, determining of three-dimensional positioning become possible in a Earth-fixed coordinate system. Therefore, by comparing the obtained positions in different epochs, displacement vector can be calculated. Displacement vectors are not suitable quantities for deformation geometric interpretation, because they must be determined in comparison with a constant basic coordinate system.
On the other hand, strain tensor can be calculated by using displacement vectors. This tensor contains invariant parameters of coordinate system and has a direct association with geometric interpretation of deformation.
Since Iran is a seismically active and has been located within the convergence zone of two rigid plates, the Arabian and Eurasian plates, It is known as a natural laboratory for studying the kinematics and dynamics of tectonic interactions. Considering these reasons, researchers aim to study the geodynamic network of Iran. The neotectonics of the Iranian plateau is complicated due to various tectonic processes, (Zagros, Alborz, Kopeh Dagh, and Talesh), subduction of the oceanic crust (Makran), and a transition zone between the Zagros fold-and-thrust belt and the subduction zone of Makran.
In this paper, using observations of GPS velocity vectors at the different station of Geodynamic network of IRAN, two-dimensional deformation of the earth crust was estimated by 2-D strain analysis within 2009 -2013. And then, in order to deformation interpretation, invariant geometric criteria like dilatation and maximum shear were evaluated.
In this regard, as the first step, coordinates of the all GPS stations was calculated in two observational epochs in global coordinate system. In addition, considering the proper map projection, the coordinates were transfered into the projection plane. In the next step, the apparent displacement vector was determined for each GPS station. Finally, using displacement vectors, strain tensor components as well as dilatation and maximum shear invariant parameters was computed for the whole Geodynamic network.
Strain tensor calculates by using displacement vector derivatives but geodetic observations, only provide the discrete values of the displacement in GPS stations. Using the data interpolation techniques such as finite element method is an appropriate way to solve this problem. In this method, the domian of the data point partitions into some smaller sub-domains and in each one, a smooth function fits in that domain by considering the specified constraints on the existing data point. In order to create a smooth surface through the data points there must be a smooth transition from one patch or element to another across all shared edges. Interpolating functions in each sub-domain must apply in some continuity conditions such as continuity of the function itself, the first or higher order derivative when passing between two adjacent elements.
Due to the nature of the Earth's crust and two-dimensional strain analysis; triangular element is the simplest geometric element in this study. In addition to determine the mathematical form of displacement function using discrete values of displacement, nonlinear finite element method was known as one of the most important point of this study. In this method, after an area meshing by Delaunay Triangular Element, C1 continuous interpolation was adopted.
The interpolating function (displacement vector) on each triangle was a rational function gained by blending three cubic polynomials that are known as the Cubic Bezier Triangular Patches. Thus, unlike the previous methods, in this case, the displacement vector in common boundary of the two adjacent Triangular elements had C1 continuity, so a smooth interpolation having C1 continuity in the entire GPS network (vertices of the triangle domain) is achievable.
Since we are working with triangles, we will utilize barycentric coordinates rather than Cartesian coordinates. So the Bezier Triangular patches (displacement vector) or Bezier triangles were defined by Bezier ordinates or control points in barycentric coordinates form.
For our interpolating problem, only numerical values of displacement vector at the three vertices of the domain triangle was provided. So the inner control points of the control mesh (unknown Bezier ordinates) must be determined. This can be done if information about the gradient or normal on the control points is given. Least Squares Minimization Method based on Finite Difference Techniques was helpful for prescribing the first order derivatives of displacement function at the vertices of each triangle for creating a C1 surface. The number of vertices around the data point is chosen by the minimum distance technique.
According to the conducted analysis, the maximum amount of shear quantity equals to 5.58×10-5 unit/year value in the southern parts of Iran. The maximum amount of dilatation quantity equals to the -2.89×10-4 unit/year value in the southeast of the country, which represents a net contraction in these areas. This contraction demonstrates the collision of the Arabian plate from the southwest to Eurasian from the Northeast, and confirms the reverse and strike-slip faults in the Iranian plateau. Also the result of the computation and the evaluations by comparison with the seismic map of the region show the success and usefulness of the presented method for deformation study of the curst.