COMPUTATIONAL MECHANICS New Trends and Applications

E. Oñate and S. R. Idelsohn (Eds.) ©CIMNE, Barcelona, Spain 1998

1

FINITE ELEMENT ANALYSIS OF THE LOAD EFFECT OF THE WATER DURING A MARKED SWELL OF THE RIO DE LA PLATA

María C. Caamaño*, Claudia L. Ravazzoli†, and Rodolfo Rodríguez‡ * CONICET. Departamento de Astrometría, Facultad de Ciencias Astronómicas y Geofísicas

(FCAyG) Universidad Nacional de la Plata

Paseo del Bosque s/n (1900) La Plata, Buenos Aires, Argentina

e-mail: crist@fcaglp.fcaglp.unlp.edu.ar, web page: http://www.fcaglp.unlp.edu.ar/ † Departamento de Geofísica Aplicada, Facultad de Ciencias Astronómicas y Geofísicas (FCAyG)

Universidad Nacional de la Plata Paseo del Bosque s/n

(1900) La Plata, Buenos Aires, Argentina e-mail: claudia@fcaglp.fcaglp.unlp.edu.ar, web page: http://www.fcaglp.unlp.edu.ar/

‡ Departamento de Ingeniería Matemática Universidad de Concepción

Casilla 4009, Concepción, Chile e-mail: rodolfo@gauss.cfm.udec.cl, web page: http://www.ing-mat.udec.cl/

Key words: Geophysics, Water Loading, Crustal Deformation, 3D Elasticity.

Abstract. This paper describes an application of the finite element method to the quantitative investigation of the deformation field of the local Earth’s surface during a marked swell of the Río de la Plata. The computational model is based on the hypothesis of linear elastic behavior for the rocks and considers the main geologic features of the upper crust. We analyze the response of the model for the normal displacement of the surface under different boundary conditions, number of layers and mesh refinements. We present numerical results showing the implementation of the method for a realistic range of elastic parameters, estimated from geophysical data.

María C. Caamaño, Claudia L. Ravazzoli, and Rodolfo Rodríguez.

2

1 INTRODUCTION As it is well known, the solid Earth is continuously deformed by a variety of processes of

external and internal nature, whose study constitutes one of the aims of the geophysical sciences. Particularly, the measurement and analysis of the deformation of the Earth due to surface traction combined with other geophysical information has been an important tool for the study of its mechanical behavior. For a long time it has been observed that large masses of water produce important depression or inclination of the physical surface even at large distances from the loads. In some cases this effect may cause temporal fluctuations in some geophysical variables such as gravity, tilt and geodetic position.

In this work we analyze the hydrostatic loading exerted on the Earth’s surface by a swell of the Río de la Plata, one of the widest in the world. Because of certain meteorological effects (such as strong southeasterly winds), this river frequently shows important increases of volume, which are observed in the registers of a continuous network of tide gauges around the coastline. That encouraged us to develop a simplified numerical model to compute the deformation field of the lithosphere in the area of the river during the forementioned swells. For observation points close to the loads (distances smaller than 30°), this Boussinesq-type problems were classically solved with convolutional integral methods based on Green’s functions for uniform halfspaces1. To overcome the physical and geometrical limitations inherent to that approach we implemented a finite element method, which allowed us to consider geological strata of variable thickness. It should be mentioned that finite element modelling was also applied in this context by some other authors such as Beaumont and Lambert2 and Sato and Harrison3.

2 CONSTRUCTION OF THE GEOPHYSICAL MODEL

Pacific Ocean

Argentina

Paraguay

Uruguay

Brazil

AtlanticOcean

Río de la Plata

Figure 1: Map showing the location of the area of study in ellipsoidal geographical coordinates.

(latitude and longitude).

María C. Caamaño, Claudia L. Ravazzoli, and Rodolfo Rodríguez.

3

To define our three-dimensional model of the lithosphere in the surroundings of the river, we first selected a two-dimensional area extending from 29° 45’ S to 39° 45’ S and from 51° 30’ W to 61° 30’ W, shown in Figure 1. Such dimensions (approximately 900 x 900 km2) were chosen so as to contain the whole surface of the river and to include points sufficiently far from it as well. This was necessary to ensure (as much as possible) the validity of the boundary conditions applied for the numerical solution of the problem.

As described in a previous work4, we analyzed the available geological and geophysical information to incorporate in our model the main subsurface geological features of the region. We considered for the lithosphere a simplified three-layer scheme, as is shown in Figure 2.

0 .0 2 5 - 7 km

3 3 km(M oh orovicic D iscon tin u ity)

1 1 5 km

Figure 2: Simplified model of the three-layer lithosphere.

The uppermost layer consists of a variety of Cretaceous-Cenozoic sedimentary fill with seismic compressional velocities lower than 5 km/sec. The thickness of this layer is highly variable and ranges from 0.025 to 7 km, as can be seen in Figure 3.

-39.75

-37.25

-34.75

-32.25

-29.75

-61.50 -59.00 -56.50 -54.00 -51.50

-7

-6

-5

-4

-3

-2

-1

0Uruguay

Argentina

Brazil

Figure 3: Map of the top of basement in ellipsoidal geographical coordinates (latitude and longitude).

The contours show depth (in km) measured from the surface of the model.

María C. Caamaño, Claudia L. Ravazzoli, and Rodolfo Rodríguez.

4

The sedimentary layer is supported by the basement, a harder rock complex mainly composed of Precambrian and Paleozoic metamorphic and Upper Jurassic volcanic rocks with compressional velocities higher than 5 km/sec. The basement layer also has variable thickness and is supposed to extend from the bottom of the sedimentary layer to a constant depth of 33 km from the surface5. The bottom layer of the model represents the lower lithosphere, for which we adopted a constant thickness of 82 km (from 33 to 115 km depths), according to Preliminary Reference Earth Model (PREM)6.

In order to take into account the effect of curvature of the Earth’s surface we worked with rectangular coordinates referred to the WGS84 system7, a widely used global geocentric reference frame. For simplicity, the real height of the topography was not included in the model, so we obtained a smooth ellipsoidal representation for the terrestrial surface in the area. This also implies an ellipsoidal shape for the Mohorovicic discontinuity (which defines the bottom of the crust) and for the lithosphere-asthenosphere contact (see the illustrative scheme in Figure 5 in the next section).

The elastic parameters of the sedimentary layer and the basement (assuming isotropy) were computed combining the velocities previously mentioned (measured by Bracaccini8 in seismic refraction surveys) with Nafe-Drake relations9 to estimate the corresponding shear wave velocities. With this procedure we established the following range for the elastic parameters λ and µ (Lamé constants) for the sedimentary layer:

λSmin = 6.15 GPa, µS

min = 0.824 GPa, λS

max = 18.33 GPa, µSmax = 21.44 GPa.

Similarly, for the basement we obtained: λB ≅ µB = 34 GPa.

However, for the numerical tests we also used for the basement λBlow ≅ µB

low = 10 GPa, according to one of the values given by Melchior10 for the crust. The mean elastic properties for the lower lithosphere were computed from PREM model, giving the following values: λL = 85.12 GPa and µL = 67.26 GPa.

To conclude the description, it should be mentioned that for the shallow water of the Río de la Plata and of the Argentine Sea, in this preliminary model we considered the same elastic properties of the sedimentary layer. This non realistic assumption avoids dealing with a fourth layer (water) in the model which, anyway should not affect significantly because of the small depth of that river and sea in the area of study.

3 COMPUTATIONAL PROCEDURE As we have previously stated, the aim of this work is to evaluate the load effect of water

during an important swell of the Río de la Plata. Disregarding any possible time dependence in the deformation process, we approached this problem by computing the static deformation of an elastic non-gravitating three-dimensional bounded domain subjected to the action of discrete forces acting on the top surface of the model (Figure 5).

So we first constructed an initial two-dimensional mesh of triangular elements with the least possible distortion (see Figure 4). The triangulation was designed so as to describe as

María C. Caamaño, Claudia L. Ravazzoli, and Rodolfo Rodríguez.

5

good as possible the shape of the coastline and the external limit of the river as well.

Figure 4: Initial two-dimensional mesh.

Next, we generated a three-dimensional mesh composed by three non-overlapping subdomains representing the geological strata described in the previous section. The initial mesh was later refined and table 1 shows the resulting number of nodes and tetrahedral elements for the three-layer case. The notation refined mesh means that each one of the initial edges was subdivided into two parts.

To analyze the influence of the bottom layer (i.e. the lower lithosphere) on the superficial displacement values we also made some numerical tests with a two-layer model which will be described in the next section.

Numbers of elements Number of nodes Initial mesh 1476 376

Refined mesh 5904 1401

Table 1. Characteristics of the different meshes used for the calculations in the three-layer case.

The differential equilibrium equations for the three components of the displacement field were solved by means of a standard Galerkin finite element procedure using linear interpolation functions.

María C. Caamaño, Claudia L. Ravazzoli, and Rodolfo Rodríguez.

6

Figure 5: Scheme of the three-dimensional computational model (in arbitrary scale) where we indicate the

different types of boundary conditions used.

To simulate the load effect of water we formulated the problem using a Neumann type boundary condition on the top given by a distribution of spatially variable surface loads taking non zero values only within the limits of the river. The hydrostatic pressure at each point was computed using real water height data measured during a pronounced swell due to meteorological causes. On the bottom boundary (i.e. the lithosphere-asthenosphere contact), we imposed the condition of zero displacement by means of a Dirichlet boundary condition. On the lateral faces we tested the algorithm with both types of condition: zero displacement and zero traction, as indicated in Figure 5.

The algebraic system of equations was solved using a conjugate gradient procedure preconditioned by an incomplete Cholesky factorization.

The computer program was written in Fortran 77 language and it was compiled and executed in a Sun Sparc Station 4. It mainly consists of the following steps:

1) mesh generation; 2) transformation of coordinates from geodetic to rectangular system; 3) finite element solution of the differential equations of elasticity; 4) identification of the displacement field at the surface and computation of its normal

component.

4 ANALYSIS OF RESULTS With the procedure outlined in the previous sections we obtained the three components ux,

uy, uz of the displacement vector field u at each point of the domain. However, for practical purposes we are interested in the displacement of the physical surface respect to the normal n to the initial undeformed state. Then, the component un at each nodal point on the surface was obtained by computing the scalar product between the unit outer normal to the WGS84 ellipsoid and the vector u.

In the following test we compare the response of the three-layer model for the initial and

María C. Caamaño, Claudia L. Ravazzoli, and Rodolfo Rodríguez.

7

the refined meshes, in both cases using a Neumann boundary condition on the lateral sides. We found significant improvements in the precision of the solution when the refined mesh was used. In particular it yields a better resolution of short wavelength features of the deformation pattern, mainly associated with geologic irregularities of the layer. This can be observed in Figure 6 (a) - (b).

-1.4

-1.2

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

(a)

-1.6

-1.4

-1.2

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

(b)

Figure 6: Contours of normal displacement at the surface in cm for λSmin, µS

min and λB, µB using a Neumann boundary condition on the lateral sides: (a) initial mesh, (b) refined mesh.

Next we tested the response of the same model, with the most refined mesh, for different

María C. Caamaño, Claudia L. Ravazzoli, and Rodolfo Rodríguez.

8

combinations of the forementioned elastic parameters. From the obtained results we observe that the area of significant deformation is located immediately under the loads and extends up to approximately 20 km from the coastline. The depression in this area is accompanied by an uplift in peripheral regions11 (see Figure 6, 7, 8 and 9). Furthermore we conclude that the maximum normal displacement ranges from -0.53 cm to -1.55 cm. This occurs at points within the Salado basin, where the sedimentary thickness and the water heights both reach their maximum values. The normal displacements are strongly related mainly to the elastic properties of the basement, and secondly to elastic properties of the sediments. This can be seen, by comparing Figure 7 (a) with 7 (b) and Figure 7 (c) with 7 (d).

-0.6

-0.4

-0.2

0.0

(a)

-1.2

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

(b)

María C. Caamaño, Claudia L. Ravazzoli, and Rodolfo Rodríguez.

9

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

(c)

-1.6

-1.4

-1.2

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

(d)

Figure 7: Contours of normal displacement at the surface in cm using a Neumann boundary condition on the lateral sides: (a) for the λS

max, µSmax and λB, µB, (b) λS

min, µSmin and λB, µB, (c) λS

max, µSmax and λB

low, µBlow, (d)

λSmin, µS

min and λBlow, µB

low.

The numerical tests performed using a two-layer model indicate that the bending of the lower lithosphere may have a strong influence on the values of the surface deformation field. This suggests that a model restricted only to the crust would not be enough to correctly simulate the deformation process in this area, as show in Figure 8 (a) - (b).

María C. Caamaño, Claudia L. Ravazzoli, and Rodolfo Rodríguez.

10

-0.4

-0.3

-0.2

-0.1

0.0

(a)

-0.6

-0.4

-0.2

0.0

(b)

Figure 8: Contours of normal displacement at the surface in cm for λSmax, µS

max and λB, µB using a Neumann boundary condition on the lateral sides: (a) two-layer model, (b) three-layer model.

Because of the large distances between the river and the limits of the model, the results obtained using Dirichlet and Neumann boundary conditions on the lateral sides do not show meaningful differences in the neighborhood of the loads (Figure 9 (a) - (b).

María C. Caamaño, Claudia L. Ravazzoli, and Rodolfo Rodríguez.

11

-1.6

-1.4

-1.2

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

(a)

-1.6

-1.4

-1.2

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

(b)

Figure 9: Contours of normal displacement at the surface in cm for λSmin, µS

min and λBlow, µB

low of the three-layer model: (a) using a Dirichlet boundary condition on the lateral sides, (b) using a Neumann boundary condition on

the lateral sides.

5 CONCLUSIONS As a conclusion we may state that the numerical model presented here gives a reasonable

approximation of the static deformation of the local Earth’s surface during an important swell of the Río de la Plata. These finite element computations show that the effect of water loading on the normal displacement at the surface depends on the tectonic structure of the surrounding area (relevant geometries and distribution of elastic constants). Differences in amplitude of up

María C. Caamaño, Claudia L. Ravazzoli, and Rodolfo Rodríguez.

12

to 40 per cent can be expected in the elastic response of various realistic, sedimentary, crust and lower lithosphere structures under the pressure of water loads exerted on extensions of a few hundred kilometers. Thus, we consider that the combination of this model with high precision gravity, tilt and geodetic position measurement may help not only to achieve a better understanding of the deformation mechanism but also to improve the estimate of the elastic properties in the area of study as well.

ACKNOWLEDGEMENTS We thank Ing. Juan Carlos Usandivaras for many interesting suggestions concerning this

work. We also wish to extend our thanks to the Departamento de Matemática Aplicada, Universidad de Santiago de Compostela (Spain) and to Lic. Claudio Brunini of the Departamento de Astrometría, Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de la Plata (Argentina), for providing us useful computer software.

REFERENCES

[1] A. Lambert, “The response of the Earth to loading by the ocean tides around Nova Scotia”, Geophys. J. Roy. Astr. Soc., 19, 449-356, (1970).

[2] C. Beaumont and A. Lambert, “Crustal structure from surface load tilts, using a finite element model”, Geophys. J. Roy. Astr. Soc., 29, 203-226, (1970).

[3] T. Sato and J. C. Harrison, “Local effects of tidal strain measurements at Esashi, Japan”, Geophys. J. Int., 102, 513-526, (1990).

[4] M. C. Caamaño, C. L. Ravazzoli and R. Rodríguez, “A numerical model for the study of the load effect of the water in the area of the River Plate”. To appear in the Proceedings of the Scientific Assembly of the International Association of Geodesy, Río de Janeiro, Brazil, (1997).

[5] L. Barrio, Modelado bi y tridimensional de campos potenciales: metodología y aplicación. Tesis Doctoral. Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de la Plata, Argentina, (1993).

[6] M. A. Dziewonski and D. L. Anderson, “Preliminary reference earth model”, Phys. Earth Planet. Interiors, 25, 297-356, (1981).

[7] Department of Defense World Geodetic System 1984, “Its definition and relationships with local geodetic systems”, Defense Mapping Agency, DMA TR 8350.2, Second Edition, (1991).

[8] O. Bracaccini, “Cuenca del Salado”. In Geología Regional Argentina, Academia Nacional de Ciencias de Córdoba, 407-417, (1992).

[9] W. J. Ludwing, J. E. Nafe and C. L. Drake, “Seismic refraction”. In A. E. Maxwell (ed.), The sea, 4, Part 1, Wiley-Interscience, New York, 53-84, (1970).

[10] P. Melchior, The tides of the planet Earth, Pergamon Press, Oxford, (1983). [11] P. Vanicek and E. Krakiwsky, Geodesy: The Concepts, North-Holland, (1986).