ICTJA PhD presentation award 2018 - Angel Valverde (Oral Presentation)

About 80 million years ago, in the late Cretaceous, the Iberian microplate was very different from the present configuration. Simplifying a little bit: the western part of past Iberia was above sea level whereas the eastern part was partially under the sea. Iberia location was between the African plate and Euroasiatic plate. Around 200 km far from the Eurasian plate, and 300 km far from the African.
In the north, there was an oceanic ridge that pushes counterclockwise to the peninsula, in what is known as the opening of the Bay of Biscay. At the same time, we had an eastward movement caused by the Atlantic ridge. As if that were not enough, the African plate and the Euroasiatic began to converge, causing a tectonic contraction in the Iberian microplate. As a result, two converging boundaries generated, one in the south, in the Betics and the other in the north, which produced the Pyrenees and the Cantabrian mountains. But why were mountains generated in the interior?

Some scientist suggest a complete bulge of the whole peninsula which generated hight topography in the middle of Iberia. With both subductions in the north and south and surface processes, we would have the actual configuration of mountain belts. Other scientist suggest that tectonic compression generated a folding of the entire lithosphere.  However, the shortening caused by the movements in the different layers of the crust can't be explained by these two theories, so it has been suggested several detachment levels that transmit efforts from the edges towards the interior, raising intraplate mountain chains.

Figure 1.

Through numerical models, we try to understand the mechanisms involved in the generation of the relief of the peninsula as well as the structure below the surface. Considering it is not possible to know exactly the past configuration of the Iberian Peninsula, we have made several models in order to explain the shortening differences of the crystal layers, as well as the overall geostructure.

We start with an initial mechanical model for a profile that goes from the bay of Biscay, through the Cantabrian range, through the basin of the Duero reaching the Central system. We also include a detachment level in the middle crust and two initial faults in the central system that come from tectonic movements before 80 million years. Afterwards, we introduce the temperature and rheologies for each zone of the model in order to reproduce the different shortening in the crust.

Figure 2.

Some studies suggest that the shortening of the upper crust is between 70-97 km in the northern part and the shortening of the lower crust of about 120 km, starting in Cretaceous times untill ~15My ago. That implies small lower crust subduction under the Cantabrian area. Moreover, the central system is estimated to have a maximum of 22 km of shortening at some points, in the last 40 My (dashed squares in Figure 5). Likewise, there is no record of shortening in the Duero basin. This fact does not imply that shortening has not been produced in the area, evidence could have been erased through time. Apart from the rheologies, a plasticity criterion is included for when there are high strain rates that generate high levels of stress where shear zones or fractures occur.

Figure 3. Topography, moho and LAB profiles for most representative models. SRTM15 is real topography of transect AB from figure 1. Grey dashed line is a previous study based in density changes from Carballo et al., that also include seismic data.

Looking into the results of the mechanical models, the one that best explains the different shortening is the M3 model, which has as a peculiarity a weak zone in the north that goes down to the upper mantle lithosphere. That generates the subduction, along with the compression tectonic movement. Model M2 has good results of shortening but the plastic criterion was easier than the M3 model and it does not have that weak zone going into the upper mantle. So there is no subduction, just thickening. Also, the structure of the material was not adequate in order to compare with the real geostructure.

Figure 4

Although shortenings coincide for M3, we can observe in figure 3 that the topography, the discontinuity of the MOHO and the LAB are very high than previous studies suggest.

When we introduce temperature gradients for the model (figure 1). The relation of density and viscosity in the model change. Now is not vertical nor linear profiles. Density changes due to the different temperature of the lithospheres that compounds the model, in this case, one is continental and the other is young oceanic. And viscosity changes due to temperature, depth, pressure and strain rate in what is call power-law creep.  Using these temperature profiles and different rheologies, strain rate accommodates in a distinct way than in mechanical models. As we can see in figure 4, those modifications change the geodynamic evolution of the model. And also shortenings are modified. But the geostructure is much more similar to the real one. 

Figure 5


Weak zones, either inherited or new, localise strain rate, increasing stress in the surroundings. Although, topograhy pattern is preserved in all models, geostructure changed quite a lot, and including temperature fix better with real one as strain rate is being localised in a different parts and in distinct way. Crustal subduction happens only when a mid crustal detachment level is predefined down to the lithosperic mantel. As we can see, plastic criteria also plays important role in the distribution of the weaknesses.

Finally, we want to say to the reader than new models with few changes in rheologies are coming soon (Figure 5 dashed black line), in order to have better results in crustal shortening. So be patient, we are getting there.

ICTJA PhD presentation award 2018 - Cristina Biete (Oral Presentation)

The influence of inherited continental margin structures on the stress and strain fields of the south-central Taiwan fold-and-thrust belt
Cristina Biete1, Björn Lund2, Dennis Brown1, Joaquina Alvarez-Marron1, Yih-Min Wu3,4,5, Hao Kuo-Chen6, Chun-Wei Ho6,7
Fold-and-thrust belts, the frontal-most part of a Mountain Range, have been studied for decades for its economic and social interest due to its importance on the exploration of natural resources. It has also been studied how and which could be the effects of incorporating in the deformation of the fold-and-thrust belts rocks that have been through a previous deformation history, however, this field is still under debate. The active Taiwan fold-and-thrust belt is incorporating in the deformation the rocks of the Eurasian continental margin (Fig. 1). This margin went through a previous deformation history that during a rifting period developed extensional faults with east-northeast orientations that now are at a high angle to the north-south grain fold-and-thrust belt. Moreover, since Taiwan fold-and thrust belt is currently active, it provides an un-parallel location to study the effects of this oblique collision on the current stresses driving the deformation and on the surface deformation.

Figure 1: Tectonic setting of the Taiwan orogen.

 In this post we want to show our research on the contemporaneous stress and strain fields in south-central Taiwan fold-and-thrust belt and how these may be influenced by the inherited structure and morphological features from the Eurasian continental margin, which are at a high angle to the Taiwan fold-and-thrust belt grain (Fig. 1).
To estimate the current stress field, we use earthquake focal mechanism, which are represented as beach balls (stereographic projection of the two possible rupture planes of the brittle deformation that produce earthquakes (see Fig. 2a for the graphic explanation of the beach balls)).

Figure 2: a) Schematic illustration of the three general tectonic regimes and the according orientations of the principle stress axes (after Anderson, 1951, and Zoback, 1992) (Modified from Barth et al., 2008 of the World Stress Map Project Guidelines). b) Earthquake focal mechanisms in south-central Taiwan study area.

From earthquake focal mechanisms recorded between 1994 and 2014 (Fig. 2b), we estimate the principal stress directions (S1, S2 and S3) and the resultant maximum compressive horizontal stress (SH), what gives a view of which forces are acting in the crust. This process is done through the inversion of earthquake focal mechanisms, from which we also obtain which are the most likely active fault planes orientations. To investigate if there are any differences in the principal stress directions in depth through the crust, we divide the focal mechanisms in three depth levels, the upper one to investigate the sedimentary cover and the fold-and-thrust belt, and two depth levels within the basement (Fig. 3). The results obtained for the stresses distribution throughout the crust are compared with the strain in the surface (Fig. 4).

Figure 3: Direction of the maximum compressive horizontal stress (SH) for each cluster at their respective depth level and colored depending on fault type.

To investigate the strain field in south-central Taiwan, we use data from the Taiwan GPS network, which is composed by the velocity vectors from the period 2005 through 2009 of each GPS stations in the study area (Fig. 4a). Since the strain is the derivative of the velocities, we obtain the strain from the velocity vectors of each station by the grid-nearest neighbor interpolation method using SSPX software. As a result, we obtain the grids for dilation, vertical rotation and shear strain rates, as well as the compressive strain orientations and the most probable shear planes throughout south-central Taiwan.

Figure 4: a) GPS horizontal velocity vectors. b) Dilatation strain rates. Blue colors representing compression and red extension. c) Summary of the stress results for the basement, the SH is shown with colored arrows, with their most likely fault planes.

The comparative between the south-central Taiwan horizontal displacement field (Fig. 4a), the maximum compressive horizontal stress azimuth (SH) (Fig. 4c) and the compressive strain orientation (Fig. 4b) show an overall similar pattern, in the north of the study area they are roughly sub-parallel to the absolute plate motion vector NW directed, whereas in the south they rotate nearly 45º counterclockwise. These rotations are produced in the center of the study area, where the Eurasian continental margin is entering in the deformation of the Taiwan fold-and-thrust belt (Fig. 4c). Also in this area the results show that the orientations of the most likely fault planes at depth and the shear planes in the surface are very similar to those orientations of the Eurasian margin faults, east-northeast striking, which are at a high angle to those faults and structure characteristic of the fold-and-thrust belt, with roughly north-south strike. The east-northeast striking active faults are typically reactivated as strike-slip faults and located in the basement, whereas newly formed faults in the fold-and-thrust belt are commonly thrusts or oblique thrusts. In the south of the study area, the results show an east-northeast oriented high shear strain area again with similar orientations to those found in the structures of the Eurasian continental margin.
Therefore, we interpret the southward change in the SH azimuth, in the compressive strain axis azimuth, and in the horizontal displacement field to be related to the reactivation of east-northeast striking faults inherited from the Eurasian continental margin.

Author’s institutions:
1) Institute of Earth Sciences, Jaume Almera, ICTJA, CSIC, Lluis Sole i Sabaris s/n, 08028 Barcelona, Spain. cbiete@ictja.csic.es. 2) Department of Earth Sciences, Uppsala University, Uppsala, Sweden. 3) Department of Geosciences, National Taiwan University, Taipei 10617, Taiwan. 4) Institute of Earth Sciences, Academia Sinica, Taipei 11529, Taiwan. 5) NTU Research Centre for Future Earth, National Taiwan University, Taipei 10617, Taiwan. 6) Department of Earth Science, National Central University, Zhongli District, Taoyuan City, Taiwan. 7) Central Weather Bureau, Taipei, Taiwan

ICTJA PhD presentation award 2018 - Mireia Peral (Poster Presentation)

Dynamics of small-scale subduction systems: a numerical and analogue approach

The theory of plate tectonics was well established around 1960. This theory describes the outer shell of the Earth as a number of thin, rigid plates (lithosphere) that are continuously in relative motion above the Earth’s mantle that behaves like a fluid. The relative plate velocities are around 5 cm/year. This implies that most of the tectonic processes and its geological consequences cannot be studied in real time since they occur on the order of millions years. With the aim of studying the dynamics of these long-term processes, analogue experiments (scaled models in a laboratory) and 2D/3D numerical simulations are commonly used. Combining computational and laboratory experiments to study a specific geodynamic process may provide a complete evolutionary study, complementing each method’s weaknesses and strengths. While laboratory experiments provide the physical realism and a high resolution, numerical models allow quantifying the physical parameters characterizing the evolution of the system that cannot be obtained from laboratory experiments.

Figure 1. Scheme of a typical subduction zone.
One of these large-scale tectonic processes is subduction. Subduction occurs along convergent plate boundaries where one plate (subducting plate) descends beneath another (overriding plate; Figure 1). Understanding subduction processes is of great interest among scientists since a large fraction of earthquakes, volcanic eruptions and ore deposits occur along these convergent boundaries. Some of these subduction zones are more complex than others involving two subducting plates in opposite direction. This type of double polarity subduction system is observed and proposed to occur in several regions of the Earth as in northern Italy or the Western Mediterranean, but the dynamics and physical parameters characterizing the evolution of such systems are poorly studied.

This work is based on 3D numerical and analogue experiments of small spatial scale subduction systems. Several single (one plate) and double (two plates) subduction models have been performed and analyzed. The objective is to compare and complement both methods to provide new insights into the analogue modelling of subduction systems and to better understand the main factors characterizing the evolution of double polarity subduction systems and related mantle flow.

Figure 2. Scheme of double polarity subduction experiment
performed in the laboratory. Plate 1 and Plate 2 subduct in opposite directions.
At large temporal scale, the dynamic evolution of the Earth’s interior can be modeled as a viscous flow problem. Accordingly, laboratory experiments consist on both linear viscous syrup and silicone putty representing the mantle and the subducting plates, respectively. For simplicity, the overriding plate is ignored. Single and double subduction models are performed in a Plexiglas tank of 150 x 150 x 50 cm3 (Figure 2). Plates are fixed at their trailing edge to enforce rollback behavior (retreating plates) and the 660 km lower mantle discontinuity is simulated by placing a fixed based at 11 cm depth. The width of the plates varies from 10 cm to 30 cm (600km to 1800 km in nature). Subduction is started by pushing down manually the leading edge of the plates into the syrup and the process continues due to the higher density of the plates with respect to the syrup.

Numerical model setup is similar to the above described laboratory experiments. Several 3D numerical simulations of single and double subduction systems with varying size of the box domain, boundary conditions, viscosity and plate thickness have been performed.

Laboratory experiments of double-polarity subduction show that trench (the contact between the subducting plate and the surface) velocities increase while trenches are approaching (phase 2) and decrease when trenches diverge (phase 3). This effect, produced by the asymmetrical pattern of the induced mantle flow, does not occur in single subduction models. Moreover, both single and double subduction models show that trench curvature increases linearly with time showing an unusual strong curvature for the wide plate models (30 cm;  ̴1800 km in nature) comparing with previous laboratory experiments of single subduction.

Figure 3. Temporal evolution (top view) of double polarity subduction model with 10 cm
wide plates carried out in the laboratory and by numerical modelling. Color arrows indicate
the mantle velocity field at 5 cm depth.
On the other hand, numerical results show that variable box sizes do not produce major differences in the evolution of a double polarity  subduction system. A box domain of 80 x 80 cm is enough to simulate accurately the laboratory experiment showing similar mantle flow pattern than in the large box model (Figure 3). Moreover, the interaction between the return mantle flow in a double subduction systems is studied by quantifying the stress and velocity field in the mantle (Figures 3&4). Our results show that two flow cells in opposite direction occur in the inter-plate region, decreasing in size during phase 2 (approaching trenches) and increasing during phase 3 (diverging trenches). Finally, numerical models of single subduction indicate that a thinner plate fits better the observations made from laboratory experiments arising the question whether the thickness of viscous plates may be modified in the laboratory during experiment preparation (Figure 5).

Figure 4. Double subduction model with 10 cm wide plates carried
out by numerical modelling. Color arrows show the mantle velocity field at
3 cm depth during plates intersection.
Figure 5. Analogue and numerical single subduction models of 30 cm wide plates at late stage of the evolution.


ICTJA PhD presentation award 2018 - Angel Valverde (Poster Presentation)

Everyone is familiar with the popular theory of plate tectonics.  It establishes the bases of the movement of the plates over a partially fluid layer called asthenosphere at around 150 km depth.  These movements generate different boundaries between the various plates. At the converging boundaries, the two plates collide and, one subducts under the other. For several reasons, the subduction is related to the composition of each plate that strikes and-or the convergence velocity.  Usually, this generates mountain chains along the subduction zone and if these plates transport continents, we have a continental-continental collision scenario.  That happened in the Himalayan area or without going far from the Iberian Peninsula, in the chain of the Pyrenees.

modified from :http://geophile.net/Lessons/PlateTectonics/PlateTectonics2_06.html

However, whether an orogen is of a certain size, or have different geological structures depends on many factors encompassed in the heterogeneity of the materials that set up the lithospheric plates during the collision.

The heterogeneity of the plates has been studied in these geological processes of a continental collision for years to better understand the characteristics of the materials that generate a type of deformation or orogen style.  Several scientist have used numerical models for these studies with different computational codes that solve the physicochemical relationships between materials. Over the years, these models have been improving thanks to the progress of the technology and the calculation capacity of computers to solve complex problems and equations Also, the level of detail or resolution for models has been improved allowing understand better the mechanism under this tectonic process. On the other hand, in laboratories scientist try to model this processes scaling different materials. These are known as analogous models. Some analogous models have been used to understand the influence of the area where both plates collide and one begins to subduct under the other.

However, we have seen that not all deformation styles that occur in numerical models and analogue models have their reflection in reality, and studies are still being carried out to clarify the factor that most influences one evolution or another. For example, some researchers argue that the folding of the entire lithosphere is a normal response to tectonic compression. Although, other researches suggest that there are levels of detachment between the different layers where stress is transmitted over long distances generating differences in shortening between the lithospheric layers.

We will focus on a model of a continental collision without taking into account the role played by the temperature (mechanical models) to try to know what is the controlling factor plus a mode of deformation or another of the lithosphere.

We consider two deformation modes; the first is known as the double-divergent orogen, where the boundaries of the mountain range are the main zones where strain rate is localised, generating structures in V called pop -UPS, which rise the ground near the suture area. The second type of deformation is the crustal folding. In this case, folds appear in both sides of the suture zone that can be very far from the subduction zone. Those folds can generate mountain chains equispaced several kilometres.

First, with a reference model and compressing the lithosphere 100 km, we verified that the initial subduction angle and changing the upper crust density were not the main factors that controlled the mode of deformation. We employ  7 different models in our study. Our results show that the main factor controlling the style of deformation is the viscosity contrast between crustal layers. This viscosity contrast can be changed directly in the lower crust (LC) or with the plastic criteria changing the cohesion.

Models M1, M2 and M7 present double vergent style of deformation. In M3 to M6 predominates crustal foldinf style of deformation. In particular, when the is a contrast higher than 100 of magnitude, the lower crust acts as a takeoff level, allowing cortical folding. The subduction of the plate is favoured reaching greater depths. And the deformation is concentrated in the crust.

But what if there was already a detachment level between the lithospheric plates that collide? We have to take into account that these processes take millions of years and the uncertainties increment when we try to look too much in the past. So, rheologies or geological structures are factors that we are not secure at all how where they in the past. 

What happens if we employ the same setup as model 2, a double vergent type of orogen, and we include two detachment levels at different depths? We will generate two different crustal folds correlated to the depth. Depending on the dimensions and position of the detachment level found in the middle crust, we have different folds near the suture zone.

If we compare the B model with the structure that the Pyrenees have, we can see some similar characteristics such as the subduction of the crust, two main faults of double V type orogen and others generated by the tectonic compression.

These detachment levels can be old faults reactivated or differences between rheology that evolve in different ways depending on the stress or strain rates deformation that affect them.


·         Folding appears when viscosity contrast between crustal layers is > 10².
·         Decoupling the crust from mantle favours crustal folding with lower crust acting as a decollement, localising the strain rate. Upper mantle in the overriding lithosphere remains undeformed.
·         Strong coupling between crust and mantle favours subduction. Strain rate localises in the upper mantle of the overriding lithosphere favouring upper mantle folding.
·         Plasticity is the mechanism that most controls the crustal mode of deformation, localising the strain rate favouring the subduction process.
·         Not all the lower crust acts as decollement.  
·         Mean topography is higher in the pro lithosphere than in the retro lithosphere.
·         Upper crust folding is related to decollement depth.
·         Decollement size influences in orogen deformation. If detachment level overpasses the subduction zone a main retro shear appears, and pop-up structures are likely to happen.


ICTJA PhD presentation award 2018 - Kittiphon Boonma (Oral Presentation)

Dependence of lithospheric slab buoyancy on composition and convergence rate: insights from a thermally-coupled kinematic model
K. Boonma, A. Kumar, D. Garcia-Castellanos, I. Jimenez-Munt, M. Fernández

Imagine a boat floating in the ocean. A single little hole will let water enter the boat. Its buoyancy is now changing. When afloat, the boat has great buoyancy, but as more and more water enters, the less buoyant the boat becomes – until it eventually sunk!

Let’s keep that buoyancy concept in mind and apply it to the case of the upper (lithospheric mantle) and lower (asthenosphere) mantle.

At a convergence zone, the cooler and less dense lithospheric mantle is subducting into the hotter and denser asthenosphere (mineral physics viewpoint).

The subducted portion of the lithosphere will experience hotter temperature and greater pressure as it subducts – causing its average volumetric density to change (denser = heavier). How does the density of this subducted lithosphere change? Well, that depends on factors such as the convergence velocity or the type of lithospheric mantle, both of which will be the key variables in this study.

Figure 1 Definition of mantle delamination
This study focuses on one of the major factors which controls the delamination (peeling of the lithospheric mantle from the crust) or subduction process – the negative buoyancy force (Fbuoy) of the sinking lithosphere. Our aim is to investigate the effect of the lithospheric mantle type and the convergence velocity on the lithospheric mantle’s buoyancy force. So, we used kinematic modelling approach to model the shortening process and subduction of the lithospheric mantle.

We have 2 main group of test parameters:

(i) Types of lithospheric mantle: Continental Lithosphere: Archon (Archean cratons, >2.5 Ga, highly depleted); Proton (Proterozoic shields, 2.5-1.0 Ga, intermediately depleted); Tecton (Phanerozoic mobile belts, < 1.0 Ga, mildly depleted), Oceanic Lithosphere: 30 Ma 120 Ma (short-lived, intrinsically denser than asthenosphere)

(ii) Convergence velocities: 1, 4, 10, 20, 30, 40, and 80 mm/yr

Slow convergence rate
At 4 mm/year, Proton Fbuoy has a negative trend down to Fbuoy=-1.25e12 N/m (114 km shortening), after which the trend starts to increase (Figure 2a and 3a). Tecton has a similar Fbuoy evolution as Proton, with a minimum Fbuoy=-9.51e11 N/m at around 103 km shortening (Figure 3b). This change of trend observed in Proton and Tecton is because a low convergence rate allows time for the slab to thermally re-equilibrate with the surrounding asthenosphere which increases the diffusion rate and, therefore, causing the slab to become more buoyant.

Figure 2 Example of effect of convergence velocity on Proton
lithosphere at: (a) 4 mm/yr and (b) 80 mm/yr.
Fast convergence rate
At 70-80 mm/year, Proton Fbuoy evolution no longer has a minimum point but instead continuously decreasing (Figure 2b and 3a). Tecton Fbuoy also continuously decrease but at a lower rate than Proton (Figure 3b). The fast convergence rate prevents the down-going slab from thermal re-equilibration (low diffusion rate), so the colder material will get pushed further down into the asthenosphere, compared to the cases at lower convergence rate, and therefore maintaining the overall decreasing in total Fbuoy (increasing in slab-pull force).

We concluded that:
1.     Thick and low averaged density lithospheric mantle: not affected by the convergence rate and always stay buoyant, implying that such mantle type does not favour lithosphere delamination or subduction.
2.     Proton and Tecton exhibit negative buoyancy force which show a tendency to, possibly, initiate delamination.
3.     Both oceanic lithospheres always subduct,
4.     Convergence velocity of ≥ 60mm/yr is fast enough to hinder thermal re-equilibration and lead to negative buoyancy – in the case of Proton and Tecton.
5.     Density contrast of < 50 kg/m3 across the LAB increase the slab’s ability to attain negative buoyancy.  

Figure 3 . (a)-(d) shows the effect of convergence velocity on the total buoyancy force. (a)
Proton (b) Tecton (c) Archon (d) Oceanic – 30 Ma oceanic (dash) and 120 Ma oceanic (solid)
lithosphere. (e) Effect of the convergence rate on the maximum slab-pull force.