Repositorio Institucional de Documentos

2013 63 María José Ramón Ortiga Flexural unfolding of complex geometries in fold and thrust belts using paleomagnetic vectors Departamento Director/es Ciencias de la Tierra Pueyo Morer, Emilio L. Pocoví Juan, Andrés Briz Velasco, José Luis Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es María José Ramón Ortiga FLEXURAL UNFOLDING OF COMPLEX GEOMETRIES IN FOLD AND THRUST BELTS USING PALEOMAGNETIC VECTORS Director/es Ciencias de la Tierra Pueyo Morer, Emilio L. Pocoví Juan, Andrés Briz Velasco, José Luis Tesis Doctoral Autor 2013 Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA MINISTERIO DE ECONOMÍA Y COMPETITIVIDAD Instituto Geoló g ico y Minero de España Zaragoza a 13 de Mayo de 2013 Emilio L. Pueyo Morer, Científico Titular del Instituto Geológico y Minero de España en la Unidad de Zaragoza y en su calidad de director de la Tesis Doctoral de María José Ramón Ortiga, certifica que la Memoria de dicha Tesis titulada “Flexural unfolding of complex geometries in fold and thrust belts using paleomagnetic vectors” se ajusta a los objetivos propuestos en el Proyecto de Tesis Doctoral aprobado por el Departamento de Ciencias de la Tierra de la Universidad de Zaragoza (18-04-2013). Atentamente Fdo. Emilio L. Pueyo Morer Instituto Geológico y Minero de España-Unidad de Zaragoza. c/ Manuel Lasala 44, 9º, 50006 Zaragoza +34 976 55 51 53, [email protected], http://www.igme.es/internet/zaragoza/ Zaragoza a 13 de Mayo de 2013 José Luiz Briz Morer, Profesor Titular del Departamento de Informática e Ingeniería de Sistemas de la Universidad de Zaragoza y en su calidad de codirector de la Tesis Doctoral de María José Ramón Ortiga, certifica que la Memoria de dicha Tesis titulada “Flexural unfolding of complex geometries in fold and thrust belts using paleomagnetic vectors” se ajusta a los objetivos propuestos en el Proyecto de Tesis Doctoral aprobado por el Departamento de Ciencias de la Tierra de la Universidad de Zaragoza (18-04-2013). Atentamente Fdo. José Luis Briz Morer Departamento de Informática e Ingeniería de Sistemas, Universidad de Zaragoza Pedro Cerbuna 12, 50009 Zaragoza, [email protected] http://diis.unizar.es Zaragoza a 13 de Mayo de 2013 Andrés Pocoví Juan, Profesor Titular del Departamento de Ciencias de la Tierra de la Universidad de Zaragoza y en su calidad de codirector de la Tesis Doctoral de María José Ramón Ortiga, certifica que la Memoria de dicha Tesis titulada “Flexural unfolding of complex geometries in fold and thrust belts using paleomagnetic vectors” se ajusta a los objetivos propuestos en el Proyecto de Tesis Doctoral aprobado por el Departamento de Ciencias de la Tierra de la Universidad de Zaragoza (18-04-2013). Atentamente Fdo. Andrés Pocoví Juan Departamento de Ciencias de la Tierra, Universidad de Zaragoza Pedro Cerbuna 12, 50009 Zaragoza, [email protected], http://www.igme.es/internet/zaragoza/ Resumen 1 Introducción La reconstrucción geológica del subsuelo es de vital importancia en muchos ámbitos como la exploración y extracción de recursos naturales (minería, hidrocarburos, recursos hídricos…) y por tanto tiene fuertes implicaciones socio-económicas. La reconstrucción en 3D consiste en la integración de todo tipo de datos geológicos provenientes de la exploración geofísica (fundamentalmente sísmica), de datos de sondeos (litología y diagrafías) y por supuesto de datos de superficie como la cartografía (Groshong, 1999). Los datos son normalmente escasos y provienen de fuentes heterogéneas por lo que es muy importante validar la reconstrucción. Los métodos de restitución son una herramienta de gran utilidad en la búsqueda de coherencia de una reconstrucción geológica. Restituir significa pasar del estado deformado al estado no deformado. Los métodos de restitución se basan en un conjunto de reglas geométricas, cinemáticas y/o mecánicas que se apoyan en una serie de supuestos sobre el proceso de deformación. El principio básico de restitución en 3D es la conservación de volúmenes durante el plegamiento (Goguel, 1952), así como la horizontalidad de las capas en el estado no deformado. Con ello se busca que los horizontes estratigráficos reconstruidos tengan consistencia geológica con un estado inicial y un proceso de deformación concreto. Así mismo, los procesos de restitución pueden ser de ayuda a la hora de predecir la deformación sufrida por una estructura. Existe un gran abanico de métodos de restitución que se aplican en 2D y 3D y que se desarrollaron inicialmente para cortes geológicos “compensados” (Dahlstrom, 1969; Hossack, 1979; etc.), evolucionaron a superficies plegadas y falladas (Gratier et al., 1991; Rouby et al., 2000; Massot, 2002; etc.) y actualmente se aplican ya a volúmenes (Moretti, 2005; Maerten and Maerten, 2006; etc.). Sin embargo, vemos que la mayoría de ellos presentan importantes limitaciones a la hora de restituir estructuras plegadas complejas (no cilíndricas y/o no coaxiales) y/o afectadas por gradientes laterales de acortamiento (rotaciones). Son estructuras que podemos encontrar habitualmente en sistemas de pliegues y cabalgamientos y que han sufrido movimiento fuera de plano (asumido como nulo en los métodos 2D) durante el proceso de plegamiento. 1 2 Objetivos El objetivo fundamental de esta tesis es añadir el paleomagnetismo a métodos de restitución existentes para mejorar su eficacia. El paleomagnetismo es el estudio del campo magnético terrestre del pasado a partir del estudio de la magnetización registrada en las rocas. Es la única técnica que permite describir movimientos respecto a un sistema de referencia absoluto, externo y global: el campo magnético terrestre. El paleomagnetismo ha sido utilizado con gran éxito desde los años 60 para la caracterización de rotaciones en zonas cinturones de cabalgamientos (Allerton, 1998; Sussman et al., 2012; etc.) y por tanto puede ser utilizado en los métodos de restitución para determinar el valor de la rotación y reducir así el número de incertidumbres. El uso del paleomagnetismo en técnicas de restitución fue recomendado a principios de los años 90 (McCaig and McClelland, 1992) y sin embargo sólo se ha utilizado hasta la fecha como criterio de corrección (Bonhommet et al., 1981; Bourgeois et al., 1997; Arriagada, 2004; Pueyo, 2000). Por tanto, la propuesta es utilizar en métodos de restitución el paleomagnetismo junto con el plano estratigráfico, puesto que ambos son una referencia conocida en los estados inicial y plegado. El plano estratigráfico determina la rotación de eje horizontal mientras que el paleomagnetismo determina la rotación de eje vertical, aportando así información complementaria. Por supuesto, los datos paleomagnéticos incorporados en el método de restitución deberán cumplir una serie de criterios de fiabilidad a los que también prestamos atención (Van der Voo, 1990; Pueyo, 2010). Concretamente, vamos a incorporar la información paleomagnética a dos métodos de restitución de superficie válidos para superficies globalmente desarrollables. No se trata propiamente de restitución 3D, pero creemos que desplegar superficies (horizontes estratigráficos) de forma correcta es el mejor punto de partida para una certera reconstrucción volumétrica. Superficies desarrollables son aquellas con curvatura gaussiana nula en todos sus puntos, superficies que han sido plegadas de forma isométrica (conservación de ángulos y longitudes) y por tanto podemos desplegarlas sin que sea necesario deformarlas (Lisle, 1992); son las que llamamos superficies flexurales. Del mismo modo, entendemos por superficies globalmente desarrollables aquellas que cumplen este supuesto de forma global a pesar de que hayan podido sufrir deformación en algunos puntos. 2 El primer método de restitución está basado en la triangulación de la superficie y parte del método desarrollado por Gratier et al. (1991). La superficie plegada está discretizada por una malla de triángulos en los que incorporamos el vector paleomagnético. Cada uno de los triángulos es abatido a la horizontal para después ser trasladado y rotado de forma que se minimice la distancia entre los vértices comunes. El dato paleomagnético nos determinará el valor de la rotación disminuyendo así el número de variables. El segundo método está basado en la parametrización de la superficie y parte del método desarrollado para gOcad por Massot (2002). La idea de la representación paramétrica es poder proyectar una superficie definida en 3D (estado plegado o deformado) a un plano en 2D (estado inicial o no deformado). Hay infinidad de posibilidades, pero por simplicidad y puesto que suponemos un plegamiento isométrico, se eligen unas coordenadas para el estado no deformado que sean rectilíneas y ortonormales. Este método está muy condicionado por la solución inicial y es aquí donde el paleomagnetismo viene en nuestra ayuda porque nos permite determinarla. Finalmente, extenderemos los métodos desarrollados a la restitución cartográfica en planta, válido tanto a escala regional cómo tectónica. No se trata de “desplegar” superficies sino de “des-rotar” la cartografía según el valor de las rotaciones verticales (VARs) obtenidas a partir de datos paleomagnéticos. 3 Metodología Para evaluar la bondad de los métodos desarrollados nos vamos a apoyar en modelos a escala. Los análogos han sido siempre de gran utilidad para tratar de comprender y explicar comportamientos y estructuras geológicas. En nuestro caso son cruciales ya que nos permiten conocer al detalle el estado deformado y no deformado y de esta forma poder compararlos con la restitución obtenida. Los modelos análogos desarrollados están basados en estructuras complejas reales que nos parecen significativas. En concreto modelizamos un pliegue cónico basado en Santo Domingo con el cierre periclinal de San Marzal y un pliegue curvo basado en el anticlinal del Balzes, ambas estructuras localizadas en Sierras Exteriores (Pirineos). Los modelos se han construido utilizando planchas de goma EVA y se han digitalizado por 3 medio de dos técnicas: 1) la fotogrametría, únicamente válida para la reconstrucción de la superficie superior y 2) la reconstrucción de volúmenes a partir de secciones obtenidas mediante un escáner de rayos X (TAC). Con la reconstrucción de estos modelos se realizan un gran número de pruebas para evaluar los diferentes métodos de restitución y los diferentes parámetros en juego. Analizamos la importancia del uso del paleomagnetismo, la sensibilidad de los diferentes métodos al mallado de la superficie, al punto de inicio y a la orientación del vector paleomagnético. Del mismo modo, se evalúan los resultados cuando el paleomagnetismo no se conoce en todos los puntos y viene definido con un cierto grado de error, tratando así de simular un caso real. 4 Línea argumental La primera parte de la tesis centra el problema. Después de la introducción (Capítulo 1) hacemos un pequeño recorrido de todo lo que se ha hecho hasta ahora (Capítulo 2). Se describen los métodos de restitución existentes para cortes compensados, superficies y volúmenes prestando atención a los supuestos sobre el proceso de plegamiento de los que parten (Sección 2.1). A continuación nos centramos en el paleomagnetismo: qué es, cómo se adquiere y qué requisitos de fiabilidad es necesario cumplir para su correcto uso (Sección 2.2). Introducimos también el uso de los modelos análogos en aplicaciones geológicas y su digitalización mediante fotogrametría, escáner de rayos X y láser escáner (aunque este último ha sido finalmente descartado) (Sección 2.3). Por último, nos centramos en el marco geológico sobre el que hemos basado el desarrollo de los modelos análogos: el cierre periclinal de San Marzal y el anticlinal curvado del Balzes situados en Sierras Exteriores (Sección 2.4). En la segunda parte pasamos a describir el trabajo desarrollado. En el capítulo 3 explicamos la metodología de los modelos análogos y los modelos concretos desarrollados. Tras una serie de pruebas, los modelos a escala finalmente utilizados han sido construidos con planchas de goma EVA sobre las que se ha serigrafiado una cuadrícula con minio a modo de sistema de referencia (Sección 3.1.3). Con el escáner de rayos X se obtienen una serie de secciones a partir de las cuales se reconstruyen las distintas superficies y su sistema de referencia utilizando gOcad, un programa de 4 reconstrucción geológica (Sección 3.1.4). Una vez reconstruido el modelo es posible calcular los tensores de deformación gracias al sistema de referencia conocido antes y después de la deformación (Sección 3.1.5). Para el modelo concreto del Bazles hacemos un análisis completo de las potencialidades de la metodología descrita haciendo un estudio de la deformación tras el proceso de plegamiento (Sección 3.2.3). El capítulo 4 es la descripción de los métodos de restitución de superficie en los que hemos incorporado el paleomagnetismo como nueva variable de entrada. El primero (Sección 4.1) es un método iterativo que trata de encajar los triángulos abatidos mediante mínimos cuadrados. El paleomagnetismo determina el valor de la rotación individual de cada triángulo, restringiendo así el número de variables y disminuyendo la incertidumbre. El segundo método está basado en la parametrización de la superficie (Sección 4.2) y el paleomagnetismo sirve para determinar la solución inicial que condiciona severamente el resultado. Finalmente se describen los parámetros de dilatación y deformación que nos ayudarán a valorar la credibilidad del método (Sección 4.3). En el capítulo 5 y 6 se llevan a cabo una serie de simulaciones para tratar de evaluar los métodos de restitución descritos junto con los distintos parámetros en juego; la sensibilidad del método en un primer lugar y la sensibilidad del paleomagnetismo en un segundo lugar. Los resultados obtenidos con la restitución se comparan con la superficie inicial y la deformación real sufrida durante la deformación (Sección 5.1). En primer lugar, se muestran la diferencia en los resultados al desplegar las superficies con o sin paleomagnetismo (Sección 5.2). Del mismo modo se evalúa si el tipo de malla que define la superficie y su densidad (Sección 5.3), así como el punto de inicio por el que se empieza a desplegar la superficie (Sección 5.4) tienen efecto sobre el resultado final. Para completar el estudio, se muestran los resultados de una restitución multi-superficie como punto de partida para una restitución volumétrica; además se comparan los resultados con una restitución 3D real exponiendo las limitaciones adicionales que ésta presenta (Sección 5.5). Adicionalmente se lleva a cabo el estudio sobre el efecto que tienen en la restitución los datos iniciales de paleomagnetismo. Se evalúa en primer lugar el efecto de la orientación inicial del vector paleomagnético con respecto a la estructura principal (Sección 6.1). No perdemos de vista que el paleomagnetismo viene definido con un cierto grado de error y se analiza su influencia (Sección 6.2). Además, es fundamental el 5 análisis de la restitución con datos de paleomagnetismo aislados y no definidos en toda la superficie (Sección 6.3). Para poder utilizar con mayor fiabilidad los métodos de restitución en un caso real donde el paleomagnetismo es conocido únicamente en estaciones puntuales, se desarrolla un método de interpolación que se basa en la suposición de partida de que las superficies son desarrollables (Sección 6.4). Se aplica de esta forma el método a un escenario más realista (Sección 6.5). En el capítulo 7 aplicamos la idea de restitución con paleomagnetismo a una restitución cartográfica en planta. Definimos las bases de este nuevo método (Sección 7.1) y mostramos los resultados para dos casos de estudio: el anticlinal del Balzes (Sección 7.2.1) y el sistema sur-pirenaico central (Sección 7.2.2). Los resultados obtenidos se evalúan en el apartado de las conclusiones proponiendo al mismo tiempo futuras líneas de investigación. En los apéndices se trata un tema de igual importancia pero un poco más alejado de la línea principal de la tesis y que está relacionado con la fiabilidad de los datos paleomangéticos. En el apéndice 1 proponemos una técnica para el procesado de los datos paleomagnéticos basada en el cálculo de direcciones virtuales como apoyo a las herramientas tradicionales. Para ello se desarrolla el programa “Virtual Paleomagnetic Directions” (VPD). En el apéndice 2 se describen matemáticamente tres posibles fuentes de error de los datos paleomagnéticos: el solapamiento de una componente secundaria, la deformación interna por cizalla y un mal control estructural del plegamiento producido en dos etapas. 6 1 Introduction 1.1 Focus Three-dimensional reconstructions of the subsurface are an important field in Earth Sciences due to their considerable socio-economic implications as exploration of natural resources (mining, oil, water, etc.). 3D reconstructions aim at providing a plausible image of the underground, which entails the integration of discrete and heterogeneous datasets: field observation (stratigraphic contacts and orientations, fault planes, etc.), interpretative map and cross-sections, seismic sections, borehole data… (e.g. Caumon et al., 2009 and references therein). Reconstruction techniques are based on geometric/mechanic laws and are designed to tackle areas with scarce and heterogeneous data. Restoration algorithms are an important tool to validate these 3D geological reconstructions of the subsurface. Each step of reconstruction must be checked and validated with geological criteria (e.g. Groshong 2006). Restoration is the way back from the deformed to the undeformed state (retro-deformation). Undoing the deformation and achieving an initial surface with geological meaning (balanced structure) is useful to validate the reconstruction of the folded structure and the deformation processes assumed. The main postulate in most restoration methods is the horizontality of the initial layers while restoration algorithms are based in several deformation processes as flexural slip or simple shear. At regional scales and crustal levels, we can assume that deformation does not change the total rock volume, at least overall (Goguel, 1952). These and other rules based on geological criteria are applied in restoration as well as in forward modeling. We deepen in restoration techniques in next chapter but we want to emphasize here the importance of a continuous feedback between reconstruction and restoration. This becomes especially important when complex deformation processes are implied and limited data is available. In addition, restoration tools may also be useful to predict deformation patterns for well characterized structures because the knowledge of deformed and undeformed states allows calculation dilation and strain tensors. 7 However, existing restoration methods do not always succeed for complex structures like non-cylindrical, non-coaxial and/or areas undergoing vertical axis rotations (out-of-plane motions). We suggest using paleomagnetic information, which is known in both the undeformed (horizontal) and deformed state, as an additional and powerful constraint to improve restoration methods and to reduce the uncertainty of the results. The use of paleomagnetism in restoration tools was recommended in the early 90’s (McCaig and McClelland, 1992). So far, however, relatively few researchers have tried using paleomagnetic information to double-check the rotation inferred from restoration methods, and hardly ever paleomagnetism is used as primary information of these tools. 1.2 Objectives In this PhD we want to show how paleomagnetism can reduce the uncertainty in restoration tools when it is used as a constraint, particularly for structures with out-ofplane motions. The bedding plane is the basic 2D reference to relate the undeformed and deformed states, but never could be a real 3D indicator. Our proposal is the usage of paleomagnetism together with the bedding plane as references known in both states. The bedding plane determines the horizontal rotation and paleomagnetism the vertical axis rotation (Fig. 1.1). Paleomagnetic vectors are the record of the ancient magnetic field at the time of the rock formation and we assume that they behave as a passive marker during the deformation process. Its original orientation can be known in the undeformed surface, and it is represented by the paleomagnetic reference vector. If we see the deformation mechanisms (Fig. 1.2), paleomagnetism allows reducing the number of variables, since it is a passive marker that may record the internal deformation and provides us with information on vertical axis rotation. Because accurate paleomagnetic data is necessary to improve results we also work on a good data acquisition. 8 Figure 1.1: Conceptual model of 3D restoration by integrating paleomagnetic data. Figure 1.2: Deformation mechanisms. Paleomagnetism helps determining the rotation. Paleomagnetism may be incorporated in many restoration tools. Particularly we center our study in geometrical surface unfolding algorithms valid for globally developable surfaces. Developable surfaces are those with Gaussian curvature equal to zero everywhere (Lisle, 1992). These surfaces in geology are stratigraphic horizons folded under flexural conditions that have minimum internal deformation. That implies surfaces isometrically folded with preservation of lengths and angles and consequently with preservation of area. By globally we mean that these constraints are valid almost everywhere but there are areas where internal deformation is possible. We can find this 9 kind of structures in the fold and thrust belts (FATs) of competent layers at upper crustal levels. In order to test the restoration methods we develop analog models of complex structures. Laboratory-scale models are based on non-coaxial structures of External Sierras (Pyrenees). These analogs are digitalized with photogrametry and X-Ray CT scanner techniques. In this way, models are completely characterized before and after deformation. This allows the calculus of strain ellipsoid of the folded surface and the comparison of the restored surface with the initial one. Additionally to the unfolding algorithms, we propose the usage of paleomagnetism in a map-view restoration technique. With this restoration we undo the vertical axis rotations (VARs) occurred during a narrow period of time. We do not pretend to do a rigorous restoration but help the location of deformation areas in the cartographic map. Premises Objectives Paleomagnetism (pmag) is reliable data known before and after deformation Incorporation of paleomagnetism in restoration techniques Pmag + bedding plane are complementary indicators in surface restoration The new constraint may help to locate the deformation FATs usually lead to globally developable surfaces Development of flexural unfolding algorithms (piecewise and parametric) Analog models are completely characterized before and after deformation Usage of analog models to test restoration algorithms VARs occur in a narrow period of time at tectonic scale Development of map-view restoration to unravel VARs 1.3 Outline We summarize in this section the main points of this PhD. In the background chapter we want to synthesize the relevant pervious knowledge. In the first section (2.1) 10 2.1.3 Surface restoration Here, stratigraphic horizons are represented with triangular meshed surfaces. Surface restoration is usually called 2.5D because only the top of the horizon is considered and not its thickness. They are usually extended to multi-surface restoration by assuming that the thickness of the layers is constant or varies in a controlled way. Rouby et al. (2000) break the process down into two separate steps: unfolding and unfaulting. In turn, unfolding algorithms are based on the two main deformation mechanisms extended from cross-sections: simple shear (Kerr et al., 1993; Buddin et al., 1997) and flexural slip (Gratier et al., 1991; Gratier and Guillier, 1993; Williams et al., 1997; Léger et al., 1997; Griffiths et al., 2002; Massot, 2002, Thibert et al., 2005). Simple shear assumes folds have internal deformation, whereas flexural slip assumes surfaces are globally developable (area and lengths are preserved) (Fig. 2.3). The unfaulting step involves dividing the region into blocks bounded by faults that can be solved with the techniques described in map-view restoration (Audibert, 1991; Rouby et al., 1993; Arriagada, 2004; Arriagada et al., 2008). Surface restoration can also be performed in one single step, using, for instance, finite elements (Dunbar and Cook, 2003) or parameterization of the surface (Massot, 2002). Concerning flexural unfolding, first methods developed were based on the triangulation of the surface (piecewise) while last methods are based on the parameterization of the surface. Gratier et al. (1991) and Gratier and Guillier (1993) divide the surface in zones defined by a network of rigid triangular elements that are rotated to the horizontal and then fitted with its neighbors by translation and rotation minimizing the sum of distances (Fig. 2.4A). Williams et al. (1997) follows the same process with a different fitting, they minimize the finite strain preserving the total area of each finite element (Fig. 2.4B). On the other hand, Griffiths et al. (2002) define a slip system (composed by a template surface, target surface, pin surface and unfolding plane) and the surface is systematically restored preserving the connectivity of the nodes (Fig. 2.4C). This technique is the one implemented in 3DMove software; unfortunately, it is not valid for structures that have suffered rotation during folding because a plane of movement is assumed and therefore, out-of-plane movements invalidate the method. Later on, several parametric approaches became to appear. Léger et al. (1997) define the unfolding process for a multi-surface in terms of parameterizations and solve 17 it by a least-squares method assuming: initial horizontality, bed-length and volume conservation (Fig. 2.4D). Massot (2002) uses a different parameterization with isometric constraints, curvilinear coordinates are orthogonal in the undeformed state (Fig. 2.4E). The commercial module Kine2D for gOcad is based on this technique (Moretti et al., 2006, Moretti, 2008). In any case, all these algorithms are designed for developable surfaces (flexural slip assumption) and when surfaces are unfoldable, the algorithm searches for the best solution although the result is never deterministic. Figure 2.3: Simple shear versus flexural slip in cross-section, surface and multi-surface restorations (Moretti, 2008). A) Flexural slip: lengths, thicknesses and area preserved. B) Simple shear: distances in the shear direction (di) are preserved, thickness, lengths and areas change. 18 Figure 2.4: Flexural unfolding algorithms with piecewise (A, B & C) and parametric (D & E) approximations. A) After flattening to the horizontal, triangles are fitted with rigid translation and rotations to minimize distances between common vertices (gaps and overlaps remain) (Gratier et al., 1991). B) After flattening, triangles are sewed together preserving its area and minimizing the strain (Willliams, 1997). C) Line length is preserved in a given unfolding direction. Node connectivity is preserved. (Griffiths et al., 2002) D) Folded and restored structure least-square minimization as defined by Léger (1997). E) Folded and restored surface using isometric constraints (Massot, 2002). 19 2.1.4 Volume restoration Real 3D restoration considers layers with thickness. The volume is represented with a grid or tetrahedral mesh. Lately, implicit approaches with relaxed meshed have been proposed (Durand-Riard et al., 2010). Last efforts leverage geomechanical approaches that meld the retrodeformational merits of kinematic balancing with principles of continuum mechanics (mass preservation and strain minimization; usually an elastic finite element model is used), without assumptions of plane strain and allowing heterogeneous fault interaction (Muron, 2005; Maerten and Maerten, 2006; Griffiths and Maerten, 2007; Guzofski et al., 2009). The drawback is that they are too dependant on boundary conditions apart from the computational requirements needed. There are also combined geometrical and geomechanical approaches (Moretti et al., 2005) and the usage of several techniques is also recommended (Lovely et al. 2012). 2.1.5 Paleomagnetism in restoration The use of paleomagnetism in restoration tools was recommended in the early 1990’s (McCaig and McClelland, 1992) as a way to tackle the restoration problem in 3D. So far, however, relatively few researchers have tried using paleomagnetic information to double-check the rotation inferred from restoration methods (Bonhommet et al., 1981; Bourgeois et al., 1997; Arriagada, 2004). These authors contrast the rotation data obtained with the restoration methods with real paleomagnetic datasets. Recently, Arriagada et al. (2008) have modified the map-view restoration method, developed by Audibert (1991) and Rouby et al. (1993), to include paleomagnetic data as primary information during the restoration process (Fig. 2.5). Although their approach is the first we are aware of that incorporates paleomagnetic data, it is still a 2D restoration method, as are two map-view methods that have been proposed involving paleomagnetic vectors (Millán et al., 1996; Pueyo, 2000 and Pueyo et al., 2004). These map-view methods (Fig. 2.5), which correct shortening estimated from cross-sections and calculate realistic shortening (using trigonometric calculus), have recently been applied in the Pyrenees (Oliva and Pueyo, 2007) and in the Rocky Mountains (Sussman et al., 2012). 20 Figure 2.5: Applications of paleomagnetism in restoration techniques. 1) Map-view concept for correction of shortening estimates in cross-sections (Pueyo et al., 2004). 2) Shortening errors as a function of the vertical axis rotations (Sussman et al., 2012). 3) Map-view applications. Two-dimensional restorations of the central Andes using two shortening models (A and B) by Arriagada et al. (2008). Colors indicate rotations of each block during restoration. A) Restored map for model with constrained block rotations (R0–15Awr). B) without rotation constraints (R0–45Anr). Arrow illustrates total displacement of 210 km in the center of the orocline. C) Restored map for model with constrained block rotations (R0–45Bwr). D) Blocks are allowed to freely rotate during the last stages of the restoration (R0–45Bnr). Arrow illustrates total displacement of 430 km in the center of the orocline. Total (Paleogene to present) displacement vectors for restoration R0–45Bwr. Paleomagnetism can be of great aid for structures with out-of-plane motions as it is the most reliable way to estimate vertical axis rotations (VARs). Of course, paleomagnetic data must be a proven and reliable record of the ancient magnetic field; the next chapter is devoted to this particular issue. 21 2.2 Paleomagnetism Paleomagnetism is the study of the ancient Earth magnetic field (EMF) recorded by ferromagnetic minerals present in most rocks of the crust. Igneous and sedimentary rocks acquire an initial, or primary, magnetization shortly after formation, which is usually aligned along the Earth’s field direction. This component may be unstable over time because of physical or chemical processes. Moreover, it can coexist with later, secondary, magnetizations and also register the ambient field direction. Thus, laboratory work is important to unravel the natural remanent magnetizations. The basic assumption to use paleomagnetism in many geological applications (e.g. detect absolute magnitudes of rotation) is that we can isolate the primary acquisition which is stable over time (Park, 1983) and it is a reliable record of the ancient field. Besides, the magnetic field is assumed to be generated by a Geocentric Axial Dipole (GAD) (Meert, 2009). Paleomagnetism is a good kinematic indicator to understand the processes of lateral transference of deformation. Following the pioneering works on plate tectonic reconstructions (see overview by Van der Voo, 1993), paleomagnetism has been increasingly used as a fundamental tool to assess the tectonic evolution of deformed areas all over the world because of its great and exclusive potential in quantifying vertical axis-rotations in an absolute way. In the last 50 years (Norris and Black, 1961) paleomagnetic data have been extensively applied to tackle tectonic problems at different scales in several orogenic systems (see overviews at McCaig and McClelland, 1992; Allerton, 1998; Sussman et al., 2012). In particular, paleomagnetism has been increasingly used as key quantitative information for unraveling tectonic deformation in fold and thrust belts and for defining the timing of the bending by its ability to determine the distribution and magnitude of vertical axis rotations (Elredge et al., 1985; Weil and Sussman, 2004; Yonkee and Weil, 2010). Together with classic structural geology analysis, reliable paleomagnetic vectors allow a spatial and temporal understanding of fold and thrust belts, including complex casestudies of non-cylindrical and non-coaxial structures. All orogenic regions have been or still are under study; Pyrenees (Oliva et al., 2010 and 2012) Alps-Carpathian system (e.g. Pueyo et al., 2007; Marton et al., 2011), Eastern Mediterranean (Mattei et al., 2007; Speranza et al., 2011), Zagros (Aubourg et al., 2008), Himalayas (e.g. Antolín et al., 2010), Andes (Roperch et al., 2011), Rockies (Wawrzyniec et al., 2007), 22 Appalachians (Stamatakos et al., 1996; Hnat et al., 2008 and 2009), Andes (Barke et al., 2007; Maffione et al., 2010), Rockies (Weil et al., 2010), the Cantabrian mountains (Weil, 2006), etc. For deformed areas it has been suggested that the remanent paleomagnetic vector might be treated as a strain marker, assuming that the magnetization is entirely pretectonic. These means that deformation has not taken place only by rigid-body rotations but by internal strain and, accordingly, paleomagnetic orientation may be modified during the deformation process. As a first approximation, it is assumed that the remanent vector behaves as a passive linear marker, rotation toward the direction of maximum extension (Facer, 1983; Cogné and Perround, 1985; Lowrie et al., 1986; van der Pluijm, 1987; Kodama et al., 1988; Stamatakos and Kodama, 1991). However, the experiments do not always confirm this simple passive-rotational behaviour of magnetic vector (Borradaile and Mothersill, 1989). Therefore, to ensure the credibility of paleomagnetic data we better consider only samples from undeformed areas or where rock internal strain can be ruled out (or assumed as lineal). 2.2.1 Characteristic remanent magnetization direction Calculation of a characteristic remanent magnetization (ChRM), as mentioned before, is a key step during paleomagnetic data processing. ChRMs are stable directions that can be effectively isolated from a given demagnetization procedure (thermal or alterning fields); subsequent application of stability tests is necessary to provide information on their geological significance. “Ideally, analytical methods (...) should be based on as much of the original magnetic information as possible, with minimal assumptions” (Kirschvink, 1980). Though, the best interpretation should always incorporate all available information (magnetic, geological, etc.). The most common technique used to separate different paleomagnetic components is to eye-ball select the relevant demagnetization interval observed in an orthogonal projection of demagnetization data (Zijderveld, 1967). Once the demagnetization interval has been selected for each specimen, directions are fitted by principal component analysis (PCA) (Kirschvink, 1980). PCA is a least-squares method to determine the linear and planar orientations of the data. Collinear points indicate the progressive removal of one magnetic vector and determine the paleomagnetic direction. 23 Coplanar points exist when two simultaneous components define a demagnetization circle1 (DC) (Jones et al., 1975; Halls, 1976). When a demagnetization routine completely demagnetizes a sample, the origin is the end-point of line through the demagnetization data; it is possible to calculate the resultant direction (R), while the other possibility is to calculate the difference direction (D) excluding the origin (Roy and Park, 1974; Hoffman and Day, 1978). Working with DCs, the paleomagnetic vector will be the intersection point between the demagnetization circles derived from different samples (Jones et al., 1975; Halls, 1976; Bailey and Halls, 1984). When multiple samples are analyzed, and where some samples provide clear end-point magnetizations and others give rise to DCs, some studies have focused on the problem of combining direct observations and intersections of demagnetization circles (McFadden and McElhinny, 1988). After obtaining an individual ChRM for each specimen, the Fisher (1953), Bingham (1974) or bootstrap (Tauxe et al., 1991) statistics are applied to determine the paleomagnetic mean vector of the site and its precision. Fisher’s parameters ( α and k) are standard in most paleomagnetic studies (Van der Voo, 1990). Other preliminary or auxiliary methods such as the Stacking Routine (SR) (Scheepers and Zijderveld, 1992) are more objective and automatic. In the SR approach, an individual mean is calculated from specimen vectors for each demagnetization step for a given site and builds a stacked demagnetization diagram for the site. Other approaches have been used to automatically fit the ChRM; linearity spectrum analysis (LSA) (Schmidt, 1982) seeks to objectively establish the demagnetization interval by means of the quality of the directions (linearity is related to the maximum angular deviation or MAD from PCA). On the other hand, the Line Find method (Kent et al., 1983) is a statistical analysis of linearity and planarity that takes into account measurement errors. So far, there is no software integrating all these methods and the transference of data between them is usually intricate. In Appendix I, we propose a new program developed based on the virtual directions to help finding the ChRM: Virtual Paleomagnetic Directions (VPD), which also integrates other methods. 1 A plane defines a great circle on an equal area projection 24 2.2.2 Reliability of paleomagnetic data in FAT belts Following the philosophy of the reliability criteria established by Van der Voo (1990) to evaluate the quality of paleopoles, Pueyo (2010) proposed some specific criteria of a paleomagnetic investigation focused on the characterization of vertical axis rotations (VAR) in an individual thrust sheet: 1) Rock, deformation and magnetization ages are known. 2) A minimum of five sites (ten is desirable) per thrust unit (10-15 specimens per site) is available. Site mean is characterized by α ≤ 10° (never > 15°) and k > 20 (never < 10). 95 3) There is a detailed demagnetization procedure isolating all magnetization components which should be fitted by PCA (Kirschvink, 1980) in which more than four steps should be involved in the calculation and MAD < 10° (never > 15°). 4) Field tests and error-control techniques (conglomerate, reversal, fold test and the small-circle intersection method) have to be performed to support the magnetization age. 5) Structural control is needed; fold and thrust geometry and kinematics should be known to avoid restoration errors. 6) The origin of the inclination error has to be identified among compaction, internal deformation and overlapping of directions by means of geometric techniques. 7) Rotations have to be contrasted to an appropriate reference in the undeformed foreland (absolute VAR) or in the nearest footwall (relative VAR). For the surface restoration methods proposed in this PhD we do not intend to calculate a VAR but only use the paleomagnetic vector in restoration techniques. Among the former criteria only points 1 and 5 do not need to be fulfilled: only predeformation acquisition must be ensured (specific age is not necessary) and paleomagnetic vector is used in situ (before any bedding correction). Table 2.1 summarizes possible sources of errors in the calculus of VARs that come from neglecting inherent assumptions about paleomagnetism in FAT belts (Pueyo, 2010). Note that, for the usage of paleomagnetism in surface restoration techniques, point 4 is not applicable because we use the vector before any correction. 25 Inclination errors, due to differential lithostatic load, are the most studied; corrections are proposed in (Tauxe, 2005). The other three possible sources of error with intrinsic structural (geometric) control are described in Appendix 2: overlapped paleomagnetic directions (C), rock volume deformation passively recorded by the paleomagnetic vector (D) and incorrect restoration of beds in non-coaxial structures (E). Assumption Source of error 1) For a given period of time, the EMF behaves as a geocentric axial dipole. A) Insufficient averaged out of the secular variations. 2) Natural mechanisms of magnetic field acquisition may be efficient to allow the ferromagnetic minerals for an accurate field orientation recording. B) Inclination flattening (shallowing). C) Overlapped directions. 3) The EMF memory may remain stable along the geological time. D) Internal deformation of the rock volume. 4) A paleomagnetic vector restored to the ancient reference system (paleohorizontal) allows quantifying the vertical axis rotations in this point (declination difference with the expected direction). E) Wrong bedding correction in complex areas where a reverse sequential restoration should be performed. Table 2.1: Error sources in the calculus of VARs. 2.3 Analog models In order to evaluate the restoration methods developed in this PhD we are going to use analog models. Analog models are really useful because they let us know the simulations (based on paper, cardboard, fabric or plasticine models) have long been performed by geologists to conceptually illustrate and understand complex structures at the laboratory scale. In particular, scaled analog models as sand-box experiments (Hubbert, 1937; Ramberg, 1981; McClay, 1990) have played an important role in 26 2.4.1.B The External Sierras The External Sierras (Fig. 2.9) are placed in the South Pyrenean sole thrust system, WNW-ESE trending and 100km long. They developed from Late Lutetian to Early Miocene (Puigdefàbregas, 1975; Arenas, 1993; Millán et al., 2000; Arenas et al., 2001), and caused the separation of the Jaca piggy-back basin to the North, from the main part of the Ebro foreland basin to the South (e.g. Puigdefàbregas y Soler, 1973; Ori and Friend, 1984; McElroy, 1990; Anastasio, 1992; Millán et al., 1995; Teixell y García Sansegundo 1995; Anastasio and Holl, 2001; Millán et al., 2000). The deformation of the Middle-Late Triassic to Early Miocene sediments were heavily influenced by the weak rheology of the Triassic evaporite deposits that served as a regional detachment horizon. The External Sierras display remarkable interference patterns between transverse (N-S to NW-SE) structures and the N-S trending Pyrenean folds and thrusts (e.g. Mallada 1878, Almera y Ríos, 1951). Millán et al. (1994 and 1995) suggests that the oblique structures were genetically related to the WNW-ESE thrust front, and also postulates that early in the structural evolution (Lutetian-Chattian), the External Sierras thrust system simultaneously developed to the south and west. Following the Chattian deformation, the kinematics changed significantly in the western and central segments of the South Pyrenean thrust system, as the development of the Santo Domingo Anticline (a regional scale detachment fold related to the emplacement of the Guarga basement thrust) and its associated south-directed thrust system progressively folded and/or truncated the earlier thrust structures. The evolution of the External Sierras involved a general clockwise rotation (e.g. Puigdefábregas, 1975; Burbank et al., 1987; Hogan and Burbank, 1996; Pueyo et al., 2002, 2003a, 2003b, 2004; Oliva et al., 2012a; Pueyo-Anchuela et al., 2012) that manifested itself in greater shortening towards the east (e.g. Soler, 1970; McElroy, 1990; Millán et al., 1995, 2000; Pueyo et al., 2004; Oliva and Pueyo, 2007a). 33 Figure 2.9: Geological setting of the External Sierras in the Southwestern Pyrenees displaying the location of the San Marzal Pericline and Balzes Anticline (mapping by Pueyo, 2000 integrating data by Puigdefábregas, 1975 and Millán, 1996). 34 Concerning the stratigraphy of the External Sierras, the lowest part of the wellexposed stratigraphic section is comprised of Upper Triassic evaporites (serving as the detachment horizon), marls and dolomites that are uncomformably overlain by Upper Cretaceous sandstones and limestones and Garumnian fluvial lacustrine facies. The shallow carbonate series of the Boltaña and Guara Fms. represent the Eocene platform that covers a great part of the south Pyrenean basin (Barnolas y Gil-Peña, 2001). In the western and central portions of the External Sierras, these Eocene carbonate deposits grade into the deposits of the Arguis Fm., consisting of azoic blue marls from outer ramp areas as well as shallow siliciclastic and carbonate facies from middle and inner ramp areas. The synorogenic deltaic sequences of the Belsué-Atarés Fm. span the Latest Lutetian to the Early Priabonian, and thin significantly to the west. Continental synorogenic strata, known as the Campodarbe Fm., is a 3000 to 4000 m thick fluvial sediment package which ranges in age from Late Priabonian to Stampian through the whole External Sierras, but is exclusively Oligocene in its western sector. Finally, the conglomerates, sandstones and siltstones of the Uncastillo Fm. span Late Eocene to Early Miocene. These strata border the southern edge of the External Sierras and lay uncomformably over the former lithostratigraphic units, recording the last compressive stages of deformation in the region. 2.4.2 Structural evolution of the External Sierras Chronology of deformation. The emplacement of the cover thrust system coetaneous with the Gavarnie basement thrust displays a remarkable diachronic character (Millán et al., 2000). This diachronism is very well-established all along the External Sierras and Marginal Ranges (South Pyrenean Central Unit) as attested by numerous syntectonic deposits. In the External Sierras, this time gap spans from Lutetian deformation (Balzes Anticline) to the onset of folding during Rupelian (Sto. Domingo Anticline). Conversely, the reactivation of cover structures simultaneous to the Guarga basement thrust affects the entire South Pyrenean basal thrust in a more isochronic fashion (Fig. 2.10). Regarding our examples, the Sto Domingo Anticline continued the recently initiated folding (with a faster pace; Oliva et al., 2012c) and the Balzes Anticline was passively (piggyback) relocated over a basal ramp-flat setting. 35 This situation is responsible for the present day plunging to the North observed in the northern sector of the anticline. Figure 2.10: Chronostratigraphic chart of deformation evidences deduced from syntectonic materials all along the South Pyrenean Front in relation to the main basement events (modified from Millán et al., 2000). 4D evolution of the External Sierras front. The existence of abundant paleomagnetic and magnetostratigraphic data has permitted to accurately control the rotation magnitude, age and even velocity in some sectors and allows for an integrative analysis of the 4D evolution of the thrust front. The diachronic emplacement of the basal thrust during the Gavarnie activity together with the very thin marine sedimentary thickness is responsible for the structuring of the first imbricate thrust system, the low wavelength and the high density of oblique anticlines in the External Sierras (Balzes, Nasarre, Tozal, Gabardiella, Lusera, Pico del Aguila, Bentué de Rasal, Rasal, La Peña, Fachar and Peña Ronquillo). The rotation activity of these formerly oblique structures, displaying systematic N-S present-day orientation, is normally coetaneous with the diachronic folding and thrusting events responsible for their genesis during the Gavarnie period (Fig. 2.11). This is demonstrated in those structures where synrotational sediments have been studied in full detail: Pico del Aguila (Pueyo et al., 2002; Rodríguez-Pintó et al., 2008), Boltaña (Mochales et al., 2012a), Balzes (Rodríguez-Pintó et al., 2013c), Mediano (Muñoz et al., 2013). Despite the onset of rotation is only partially established, a well36 defined diachronism has been demonstrated for the end of the rotational movement. The rotation laterally vanishes at an averaged rate of ≈ 5 km/M.a. although this velocity may be faster if a non-steady scenario is considered. Comparable values could be expected for the rotation onset as regards of the remarkable similarities found with the lateral migration of the deformation along the External Sierras front (Millán et al., 2000) or the westwards onlap of the turbiditic trough (Labaume et al., 1985). Younger rotational activity (related to the Guarga emplacement) cannot be completely ruled out but, if exists, is expected to be very small. This observation, apart from the existent paleomagnetic data, is supported by the very quick lateral expansion of the Guarga thrust front (≈ 30 mm /year), which, in turn, implies a very small lateral gradient of shortening (and equivalent associated rotations of the thrust front). Due to this complex deformation pattern (imbrication, obliquity, diachronity and two main deformation events), the External Sierras are an excellent natural laboratory to study the 3D geometry and kinematics of complex structures caused by non-coaxial axis of deformation and vertical axis rotations. Figure 2.11: 3D diagrams (not to scale) showing a schematic rotation model of the evolution of the Western and Central sectors of the External Sierras front (Pueyo et al., 2002). A) t1 (40 Ma); at the beginning of the deposit of the marl sediments (Arguis Fm.), illustrating the effect of the onset of rotation in the hanging wall of the South Pyrenean basal thrust. B) t2 (36.5 Ma); until the end of the deposit of the transitional sediments (Belsué-Atarés Fm.), which also represents the end of the rotation of the basal thrust in the Central sector. C) t3 (26 Ma); during this period, the studied area did not show any significant rotations but experimented important translation towards the south, while for the same time span, the western sector of the External Sierras suffered important rotations. 37 2.4.3 San Marzal Pericline The San Marzal Pericline is the lateral termination of the Santo Domingo Anticline, located at the westernmost sector of the External Sierras (South Pyrenean sole thrust). This deca-kilometric and apparently cylindrical anticline accommodated most of the shortening in that area. It strikes WNW-ESE, detaches along the incompetent Keuper facies and depicts parallel near-vertical limbs (Millán et al., 1995). It was active during Late Oligocene-Early Miocene (last stage of the structural evolution of South Pyrenean thrust front). Recent magnetostratigraphic studies (Oliva et al., 2012c and in prep.) in the southern flank of the anticline as well as and the reinterpretation of previous sections (Hogan 1993; Arenas et al., 2001) identify two distinct folding periods in relation to the Gavarnie and Guarga emplacements. These two periods display very contrasted folding velocities (Fig. 2.12). Figure 2.12: Kinematics of the Sto. Domingo Anticline deduced from the magnetostratigraphic studies in its southern flank (Luesia section by Oliva et al., 2012c and 2013 in prep) The Mesozoic beds as well as the marine and lowermost continental strata of Tertiary age involved in the Santo Domingo Anticline describe a cylindrical closure at the western termination of the External Sierras: the, so-called, San Marzal Pericline which folds axis orientation is 305º, 67º (Fig. 2.13). The underground western geometry of the fold reflects a quick diminishing of the plunge of the axis (Oliva, 2000; Oliva et al., 2012a). 38 Figure 2.13A: Geological setting of the Sto. Domingo Anticline (red rectangle) at the western end of the south Pyrenean sole thrust. Geological sketch map displaying the location of paleomagnetic data and cross-sections (modified from Puigdefábregas, 1975, Millán, 1996 and Pueyo, 2000). Paleomagnetic rotations by Hogan (1993-blue) and Pueyo (2000-white) are also displayed; cone axis is the mean paleomagnetic declination and its semi-apical angle represents the confidence angle (α95). Figure 2.13B: Stereographic projection of bedding poles; San Marzal Pericline and Sto Domingo Anticline. A cylindrical best-fit (Bingham’s [1974] statistics) performed with the Stereonet program characterizes the fold trend and plunge. Stereographic projections using Stereonet (Allmendinger et al., 2012 and Cardozo and Allmendinger, 2013). 39 Figure 2.13C: Balanced sections in the western termination of the External Sierras; Cross section-I: Isuerre (Oliva et al., 1996 and 2012a) and cross section-II: San Marzal and III: San Felices (Millán, 1996). Note the effect of the fold axis plunge (cone generator trend) on the geometry of the preCampodarbe sequence. 40 Figure 2.13D: Conceptual model for the Sto. Domingo anticline at the western end of the south Pyrenean sole thrust (Millán et al., 1992 and 1995). The interest of this conical structure underlies in the large rotation magnitudes associated to its genesis, as originally proposed by simple analog modeling (Millán et al., 1992). Paleomagnetic analysis, carried out in fourteen marly sites (Arguis Fm.) at both flanks of the anticline as well as in the fold termination, attests a significant clockwise rotations (CW ≈ 45°) at the northern flank and almost 20° counter clockwise rotation (CCW) at the southern one (Pueyo, 2000; Oliva et al., 2012a; Pueyo-Anchuela et al., 2012). The sites located around the fold hinge (San Marzal area) display variable and gradual rotations between both extreme terms. This particular geometry is probably due to three combined mechanisms: 1) The general lateral gradient of shortening associated to the emplacement of the External Sierras thrust system responsible for the about 30° CW rotation in average (McElroy, 1990; Millan, 1996; Pueyo et al., 1996 and 2004: Oliva and Pueyo, 2007a). 2) The lateral disappearance of the detachment level (Keuper facies) to the West, as demonstrated by the borehole records (Aoiz, Roncal and Sangüesa wells; Lanaja, 1987 and balanced cross sections by Oliva et al., 2012a) that would have produced a pinning effect and the conical geometry of the fold. 3) The southern flank CCW rotation seems to be local and it dies out to the Southeast of the Santa Engracia fault (Pueyo et al., 2003). Consequently this rotation looks to be a local effect probably caused by the pinning of the fold and the impossibility of the northern flank in accommodating more rotation. 41 Therefore the Santo Domingo anticline and the San Marzal lateral termination seems to be a pivot-point conical fold (in terms of Allerton, 1994) and it is a perfect paradigm of a complex-structure that can be used to test the reliability of our restoration method. 2.4.4 Balzes Anticline The Balzes Anticline (BA) (Figs. 2.9 & 2.14) represents the southeasternmost structure of the External Sierras in the Southern Pyrenees. It is a 17 km-long curved structure with a fold hinge trending N011E in the northern sector passing to N152E in the southernmost sector, therefore, in map view, it displays an apparent bending of about 40° (southwestwards convex). Another interesting aspect is its connection with the N-S Boltaña anticline to the N, which partially overlaps the Balzes fold axis to the East. In contrast to the western sector, the particular stratigraphy involved in the BA comprises three main marine platform sequences during the Eocene (de Federico, 1981; Barnolas and Teixell, 1994): the Ypresian Alveoline limestones of the first platform outcrop in the core of the structure; the Boltaña Formation (late Ypresian, locally Cuisian), the second platform, represented by ≈300 m of shallow limestones and siliciclastic input; and the third platform, the Guara Formation, made of up to 650 m of Lutetian limestones. The sedimentation of the Guara Formation was determined by the growing of the Balzes Anticline as attested by an angular unconformity (Fig. 2.15) observed in its western flank (Millán et al., 2000; Barnolas and Gil-Peña, 2001; Rodríguez-Pintó et al., 2012b). On top of the Lutetian, in the northern part of the anticline, the Sobrarbe deltaic formation (Bartonian) marks the transition to continental conditions, the onset of which is clearly indicated by the thick Campodarbe Group (Puigdefàbregas, 1975) cropping out in the core of the vast Guarga Syncline (Jaca Basin). 42 Figure 2.20: Balzes Anticline geometry and kinematics constraints (Rodríguez-Pintó et al., 2013c). A) Balzes curvature. Bending diagram (VAR versus structural trend) in the anticline, data from Boltaña (BA) and from Pico del Aguila (PAA) anticlines are also shown. B) Rotation velocity in the Balzes anticline. Paleomagnetic rotations derived from mean (robust) values obtained for discrete temporal gaps (2-3 Ma). 49 50 3 Analog models Laboratory scale models are designed to evaluate the goodness of the restoration method. So far, we have developed static models in which only initial and final states are detailed. With these experiments we can fully characterize the geometry of any structure before and after deformation. Thus, they are perfect tools to compare the restored surface obtained with the restoration method with the initial and expected one. The analogs we want to simulate are complex structures from fold and thrust belts (FAT belts) at upper crustal levels. In this work we have only modeled folded layers without considering the faults; the whole structure is divided by regions bounded by faults and each block is treated individually. In the upper crustal levels (within the first 5-6 km depth), competent layers such as limestones or sandstones behave more typically with a flexural slip mode (Ramsay and Huber, 1983). The preservation of lengths and angles (and consequently also areas) during folding is called in differential geometry as isometric bending. A remarkable property of isometric bending is that Gaussian curvature (the product of the two principal curvatures) is invariant and equal to zero everywhere; the surfaces are developable (Lisle, 1992). However, there may be some localized deformation in specific areas, i.e. flexural flow on the flanks or tangential longitudinal strain at the hinges (Ramsay, 1967). Therefore, we can speak about globally developable surfaces. The other important property we want to simulate is paleomagnetism. We are going to plot lines on the surfaces to represent paleomagnetic vectors (declination component). The initial paleomagnetic vectors in a real scenario are not contained in the surface as those simulated, but we can simply rotate them to have null inclination and then be embedded in the surface. Only the declination is relevant. It is worth saying that, we assumed a perfect primary record of a GAD magnetic field. Additionally, an interesting aspect that analog models may offer is the possibility to quantitatively measure the internal deformation, and not only qualitatively evaluate the results. An appropriate reference system allows the quantification of displacement, area or volume change and strain. The reference system proposed is an orthogonal grid drawn on each surface. The comparison of nodes location of two adjacent will allow us to reconstruct the strain ellipsoid (see more details later). 51 Otherwise, to digitalize or reconstruct the analog models we use two techniques already introduced: X-ray computer tomography scanner for volumes (built from a dense set of serial cross-sections) and photogrammetry for surfaces (built from a set of referenced nodes). The first technique requires materials with radiological contrast. The restoration methods developed in this work are for single surfaces and do not for the whole volume. In that sense, only reconstruction with photogrammetry would be enough. However, the technique of CT scanner allows evaluating inner surfaces and opens a wide range of possibilities for future researches, like understanding complex structures and characterizing its 3D deformation patterns, as well as validating others 3D reconstruction and restoration methods and software. We show its potentiality with the example of the Balzes Anticline. 3.1 Methodology In this section we describe the particularities of the analog models developed. These models are valid for flexural folds, they incorporate paleomagnetic data and an orthogonal reference system, and are appropriate for X-ray CT and photogrammetry reconstruction. We first describe the principles of photogrammetry because it is a technique with few requirements, and we focus later on the needs for the X-ray computer tomography (Ramón et al., 2013). The upper surface of models built for CT reconstruction can be also reconstructed using photogrammetry. In the last subsection we explain how to calculate the strain ellipsoid to quantify deformation using the orthogonal reference system. 3.1.1 Photogrammetry: principles and settings Photogrammetry is a technique that only requires a few photographs taken from different angles to reconstruct a 3D model. Any conventional camera is valid as long as all photos are taken with the same focal length. The software used for the reconstruction is PhotoModeler. This program uses the description of the camera (including data on the focal length, imaging scale, image center and lens distortion) to build a proper geometrical relationship between points on the photograph and points in 3D space. This 52 information is obtained by the camera calibration, performed by taking a minimum of six photographs of a reference grid provided by de program. Once the camera is calibrated we take photographs of the model from different angles, in order to see all points from at least two views, and recommended from at least three. Then, we proceed with the referencing process, marking on two or more different photographs points that represent the same physical object in space. We mark all points that define the grid (the reference system that we have drawn on our model). From these set of xyz points we build the triangular mesh that defines the surface. Accuracy depends heavily on the precise marking (or not) of locations on the images. A normal relative accuracy of 1:5,000 means that for an object with a 1 m largest dimension, PhotoModeler can produce 3D coordinates with 0.2 mm accuracy at 68% (one standard deviation) probability. Figure 3.1: Analog model reconstruction with photogrammetry. Reference points are marked from several photos to build the 3D model (right). All photos must be referenced between them and all points must be seen at least in two photos taken from different angles. 3.1.2 X-ray CT: principles and settings Conventional CT scanners used in medicine usually have millimeter-scale resolution based on using low-energy X-rays (below 125 kV). A scanner with these characteristics is sufficiently powerful for our purpose. If the object scanned has a low 53 radiological density, it is possible to increase the intensity of the emission source. This should, however, be avoided when possible because of the risk of gantry overheating (in particular, the X-ray tube). Apart from the reological considerations, the proper selection of the materials to build the model is an important factor in achieving clear CT images. We used a General Electric HiSpeed FX/i CT scanner at the Royo Villanova Hospital in Zaragoza (Aragon Health Service, SALUD) (Fig. 3.2) in collaboration with L.H. Ros and the Radiology Service technicians. For the current study, the settings selected to optimize the digital reconstruction were (Fig. 3.2): 1) axial scans, rather than helical, because they generate a sharper image, with the minimum slice thickness allowed (1 mm); 2) slices spaced 0.5 to 1 cm apart, which is close enough for the required resolution; 3) a high resolution chest CT protocol with a beam energy of 120 kV and low current of 180 mA to avoid gantry overheating; 4) the lung window to properly view the images on the CT system; and 5) the DICOM format to export the data to the 3D reconstruction software. Figure 3.2: General Electric HiSpeed FX/I CT scanner at Royo Villanova Hospital in Zaragoza. Left bottom inlet. Scanner settings in the main menu. 54 3.1.3 Materials and rheology Many of textile fabrics accomplish the basic geometric principles of isometric folding. They can be complexly folded but they will always be (globally) developable. Therefore, simple fabric simulations can be a very useful basis for realistic 3D models (Fig.3.3). Among the many possible materials, ethylene vinyl acetate sheets (EVA, also known as expanded rubber or foam rubber) allow almost any kind of deformation, including flexural flow and flexural slip as well as tangential-longitudinal strain. It has an amorphous structure, its density can vary widely, from 50 to 200 kg/m3, and it has very low water absorption (≈0.07%). On the other hand, its tensile strength ranges from 2 to 10 N/cm2, while tear strength varies between 2 and 4 N/cm2, and it can be elongated by as much as 500 %, although common values are around 200-300 %. Due to its versatility in industrial applications, it is commercially available in thicknesses covering three orders of magnitude (0.5 to 500 mm). The stacking of a given number of EVA sheets (stuck together or free to move) also allows the stratigraphic thickness of the model and the expected mechanism of deformation to be varied. Unlike other materials, EVA’s radiological contrast is high enough, and its boundaries can be imaged with sufficient clarity to be accurately redrawn by the image processing software. Alternatively, ethylene propylene diene monomer rubber (EPDM [M-class] rubber) can also be used. Alternation of different mechanical properties, thicknesses and cohesion between layers produces infinite possibilities and allows modeling any case under study. The thicknesses of the modeled stratigraphic pile and the wavelength of the folds have to be adapted to the limited size of the CT scanner (usually less than 60 cm in diameter) and the circular geometry of the CT sections has to be taken into account in the model design. 55 Figure 3.3: Modelization and scanner of several analogs. 56 In addition to the layer modeling, the reference system characterization is of crucial importance. This reference system is defined with two sets of orthogonal parallel grid lines. We tested many different kinds of materials to simulate the sets of lines: plastic and metallic meshes, strings and cords, and even cloth with linear relief. In most cases, however, there was insufficient radiological contrast, while in the case of metallic lines the absorption was so high that the radiological image was over-exposed. Moreover, the use of specific layers such as plastic meshes to simulate the set of grid lines has mechanical consequences: 1) it implies detachment from the underlying bed; 2) it does not allow the effect of deformation on the grid lines to be quantified; and 3) the change of the rheology has also consequences for the final geometry. Solid linear elements (cords, strings, wires, etc.) cannot effectively be stuck to an EVA sheet (different mechanical properties) and if they are free to move, they cannot be trusted by definition to provide an accurate reference system. Similarly, the use of cuts or marks in EVA sheets also affects the rheology. Therefore, to define the reference system we decided to paint a set of parallel lines on the sheet using highly X-ray absorbing inks, liquids and paints. A wide variety of materials were tested for this purpose, all of them having high electrical conductivity (Fig. 3.4A), including graphite, gold and silver inks, aluminum paints, various types of glitter, etc. Of these, minium (red tetraoxide lead) paint was, by far, the most successful material. We found it gives a sharp radiological signal without serious streak artifacts in the CT image. Other potentially suitable materials such as graphite, gold, silver and aluminum had too little mass to absorb enough X-ray photons and hence were not detected in the signal or they were but only in a very faintly way. The minium was screen-printed onto the EVA sheets, which can be done with a good accuracy. Indeed, current screen printing technologies allow computer-aided design of the screen mesh (usually made of nylon and polyester) and mesh sizes up to 0.5 mm. Minium can be used in place of the usual screen printing inks without affecting the performance of the printing process. The grid that constitutes the reference system (needed to monitor internal deformation in 3D) is made by two orthogonal sets of parallel lines. We explored screen printing the second set of parallel grid lines using a mixture with a different ratio of minium and turpentine or, alternatively, a different quantity of ink. The amount of ink can be varied by changing the line width and an advantage of this approach is that the 57 minium/turpentine mixture for the main and secondary sets of parallel grid lines can be applied at the same time. We found that varying the line width between 0.5 and 2 mm in the screen worked well. These different amounts of minium produced sufficiently different intensities to allow identification of each set of lines in the CT images (Fig. 3.4B). Figure 3.4: A) CT image of lines applied with the different conducting materials tested. B) CT image of lines containing different proportions of minium (left to right) to reproduce a secondary reference grid. Note that the signal from the minium mixture does not change from the top to the bottom section. C) Screen-painted EVA plates using dissolved minium. Different “ink” width (dash and continuous) trials to test the different radiological brightness. 58 the vertices of the tetrahedron. There are three tetrahedra for each triangular prism (Fig. 3.9) and we calculate its mean ellipsoid in order to have a mean value for the whole volume of this triangular prism. Figure 3.9: Volume tessellation Dilation is the change of volume between initial and final tetrahedron: ( ) ( ) ( ) 1detdet − ′ =−= VVvolumevolumevolumedilation initialinitialfinal Using basic properties of the determinant, it can easily be shown that )det(/)det()det( VVM ′ =. Moreover, ( ) 2 )det()det( − =MAellipsoid , and thus () 11 1 )det(1)det( 321 321 2/1 −⋅⋅=− ⋅⋅ ==−= −kkkAMdilation ellipsoid λλλ To describe the anisotropy of the ellipsoid, we calculate the P’ and T parameters defined by Jelinek (1981). We are interested in the ratios of the semi-axes, rather than in their differences. As a measure of their scatter, we consider the normalized anisotropy: ()()( [ ) ] 2 3 2 2 2 1 2exp' ηηηηηη −+−+−⋅=P where ηi are the logarithms of the semiaxes )ln( ii k= η , and η their mean value 3)( 321 η η η η + + = . The shape factor, defined as 31 312 2 ηη ηηη − −− =T, characterizes the shape of the ellipsoid. An ellipsoid is said to be rotational prolate (prolate, neutral, oblate, rotational oblate) when ( ) 121222323 , ~ , ~ , ~kkkkkkkkkkkk =<<=<<= , and thus T=-1(-1<T<0, T=0, 0<T<1, T=1). 65 3.2 Analog models With the analog models we do not pretend to accurately reproduce the structure of San Marzal and the Balzes Anticline, neither to do an exhaustive analysis about them. These structures are selected because of its complexity (a conical fold and an arched anticline with a related conical fold in its inner arch). Their geometry is a good case to test the capabilities of the restoration methods and the CT scanning. In any case, and following the philosophy of analog modeling, our models assume an evaporitic core that has been modeled by air in the core of the anticline (a reasonable assumption considering the rheology). The cover rocks on top of the Middle Eocene platform (Boltaña Fm) has not been modeled, according to the deformation ages (synfolding). The model was designed as a “static” reproduction of the evolution of the anticline (we apply the finite deformation at once, without considering the actual kinematic). Several models have been developed based on this two complex structures (Marzal and Balzes) using different materials and scales (Fig. 3.10). It has been a laborious work in which we have learned the modelization and the reconstruction technique. At the end, we have selected two of the best reconstructed analogs to show the results of this technique and the restoration methods. We have reconstructed the San Marzal model using only photogrammetry and the Balzes Anticline with both techniques. Figure 3.10: Analog modeling of several models of San Marzal and Balzes. 66 3.2.1 San Marzal Pericline The San Marzal closure of the large-scale Sto. Domingo detachment anticline is modeled considering the pre-Campodarbe sequence, and particularly, the Guara Formation. The air within the anticline core represents the pre-Guara formations including the Triassic evaporites. The analog reproduces five main structural features (Fig. 3.12): 1) approximately the San Marzal relation between the wavelength (≈ 8 km) and thickness of the modelled layer (≈ 180 m); 2) strong immersion of the fold axis; 3) pseudo-parallel flanks in the Sto. Domingo Anticline; 4) an approximately 45° clockwise rotation in the northern flank; 5) the southern flank is assumed to be in structural continuity with the Ebro foreland basin and it will be used as the pin-line. Figure 3.12: Analog model of the San Marzal Pericline. A) San Marzal cartography (from Fig. 2.13A). Approximate wavelength and thickness of the modelled layer displayed. B) Orthophoto of the analog model. Approximate wavelength and rotation of the northern flank displayed. C) EVA foam model on which lines parallel and perpendicular to the paleomagnetic reference vector have been screen-printed. D) Final 3D reconstruction of the model with the z coordinate displayed. 67 The model is built with an EVA plate of 0.3 cm thickness. Therefore, the relation between the wavelength and thickness is approximately 50 (15 cm / 0.3 cm) and is similar to one of the real structure ≈ 44 (8 km / 0.18 km). The paleomagnetic vectors (projection of the declination component) are featured as lines printed on the EVA surface. A rectangular grid screen-printed on the surface provides two possible references. The 2020 × grid is composed of rectangles of 25.1 × cm (Fig. 3.12). We have reconstructed the San Marzal model using only photogrammetry because after many attempts we have found the CT reconstruction impossible. For a plausible CT reconstruction, lineations must be as perpendicular as possible to the cross-sections. Because of the strong folding of the model (overturned limbs) and the strong differential rotation between the flanks we were unable to fulfill this premise. 3.2.2 Balzes Anticline We model the upper and better exposed thrust sheet involved in the Boltaña Formation (unit 2 in Fig. 2.19). The evaporitic core is modeled by air and the syn-Guara Formation has not been modeled. Thus, the top surface of the model represents the base of the Guara Formation and the bottom surface the base of the Ypresian limestones. The geometric scaling of the model obeys some key features: 1) the real variation in the fold axis trend (stereoplots in Fig. 2.18), 2) the relation between the wavelength of the anticline (≈ 6 km) and the thickness of the modeled stratigraphic pile (≈ 300 m), and 3) the differential vertical axis rotation between the northern and southern sectors (27°). Generating an oblique structure based on the Balzes Anticline, we see how a secondary fold is formed in the inner part, which could correspond to the Boltaña Anticline southern termination near Paules de Sarsa. Using all this information the analog model is built with two EVA sheets of 58 x 38 x 0.5 cm glued together, giving a total thickness of 1 cm. In this case the relation between the wavelength and thickness is approximately 16 comparing with the approximately 20 of the real one. The EVA sheets are screen-printed with a squared grid of 1 x 1 cm and line widths of 1 and 1.5 mm. One sheet is screen-printed on only one side and the other on both sides, thus we model three surfaces that represent the base, the top and the middle (and neutral) surface of the stratigraphic pile under study. 68 First of all, we scan the EVA sheets in an undeformed state (horizontal) to set up the reference system. Subsequently, we deform the sheets following the Balzes kinematic model and scan it again (Fig. 3.13). We then reconstruct the model from the DICOM cross-sections in gOcad. We also check the reconstruction from CT images against the reconstruction of the upper surface obtained using photogrammetry. Figure 3.13: Analog model of the Balzes Anticline. A) Balzes Anticline cartography (Fig. 2.19). Approximate wavelength displayed. B) Orthophoto of the analog model. Approximate wavelength displayed. C) CT scanning of the analog model. D) Reconstruction with gOcad of the analog model from the CT images. 69 3.2.3 Analysis of Balzes Anticline CT model With this example we want to show that these models can be very useful for understanding certain features related to complex settings: A) folded lineations, B) the 2D distribution of deformation ellipses on different surfaces within the model, and C) the 3D distribution of strain tensors in the folded volume. 3.2.3.A Lineation analysis This technique is valid for studying any passive geological lineation. In this section, we compare the grid lines of the model with paleomagnetism, but we could study any other linear element (paleo-current, stress, etc.). Paleomagnetic vectors (their local projection on the model surface; that is to say, only the declination information) can be seen as a type of lineation (Sellés, 1988; Stewart, 1995; Pueyo et al., 2003) and are feasible structural markers that can be clearly established for both the pre-and postdeformation states. In this case, paleomagnetic data are assumed to be pre-folding. First, we need to rotate the entire model to converge to the real axes of the structure. As the reference set of lines in the model had no inclination, we need to apply a rigid-body rotation to our contrived paleomagnetic record (one of the sets of lines) to converge with the real dataset. Inclination has been modeled to fit the real data (≈40°) rather than the Eocene reference expected in the Pyrenees (53°). Now, we consider lineation patterns separately in the northern and in the southern flank of the anticline. We select specific sites on the model simulating an outcrop on each side of the anticline. These data are approximately at the same structural location as the real dataset. We project paleomagnetic vectors before and after bedding correction (Fig. 3.14A). Paleomagnetic data before any correction are clustered into two groups corresponding to the western and eastern limbs of the anticline. As expected, data is grouped after the bedding correction (ABC), and the clustering is better than seen with the real noisy paleomagnetic data (Fig. 3.14B). Given the secondary origin of the fold curvature (that we applied to the model), there is a 22° difference between the mean paleomagnetic direction in the northern and southern sectors. This difference is slightly smaller (4°) than detected in the real dataset. These small errors in the lineation, as well as those highlighted by the bedding poles, are not unreasonable. They are likely caused by the analog model, which does not 70 perfectly reproduce the natural geometry. Once again, the fold axes of the sectors are similar but not exactly equal to those calculated for the real structure. Figure 3.14: Paleomagnetic analysis of northern and southern flank. Stereographic projections of paleomagnetic data before and after bedding correction (BAC, ABC) as well as the bedding plane (S0) for each site with the calculation of the fold axis: A) Analog model; and B) Real data. 71 3.2.3.B Surface analysis Before considering the internal deformation, we analyze the results of the uppermost surface, comparing the results obtained from the reconstruction of the CT images with those from photogrammetry, the complementary technique we use for surface reconstruction. Our aim with this was to validate the 3D reconstruction based on the CT images. First of all, we compare the distribution of dilation (=[Areafolded– Areainitial]/Areainitial) across the entire surface. The second technique is more accurate for the upper surface reconstruction and gives clearer results but they are equivalent in meaning (Fig. 3.15A). This is basically due to the greater accuracy of the photogrammetric method in reconstructing the exact location of the nodes (intersections between the two sets of lines). If we now analyze all the three surfaces derived from the CT modeling (Fig. 3.15B), we observe tangential-longitudinal strain in line with Ramsay (1977) and Gairola (1978), extension in the anticlines outer hinges of the upper surface (positive dilation), conservation of area in the middle neutral surface and compression (negative dilation) of the inner hinges in the lower one. We observe consistent senses of dilation in the synclines: positive on the outer arc and negative on the inner one. It is worth noticing that the areas with clear dilation are significantly different at the two surfaces (upper and lower): in the lower surface the maximum dilation is concentrated in a much smaller area and is more intense than in the upper one. This observation fully agrees with the expected differences in arc lengths between outer and inner hinges. The most deformed areas are precisely the ones of maximum mean curvature (mean between the two principal normal curvatures: M=(k1+k2)/2). Since we are unable to derive the anisotropy from the curvature, we plot the strain ellipses in terms of the ratio of the axes and the orientation of the main axis. In the CT model, the areas of higher anisotropy (ratio between major and minor axis) correspond with areas of compression: the syncline between Balzes and Boltaña anticline in the upper surface (synclastic synform according to Lisle and Toimil, 2007) and both anticlines in the lower surface (synclastic antiform). This observation agrees with the distribution of dilation across the upper and lower surfaces. As suggested above, this may be caused by the concentration of deformation in the inner-arc zones, which have to accommodate an equal amount of volume change in a smaller deformed volume. On the other hand and, as it would be expectable, the Gaussian curvature (G=k1·k2; Gauss, 1827) is concentrated in the area 72 of superposed folding (Boltaña-Balzes) and there are no significant differences between the three surfaces. 73 Figure 3.15 (continued): A) Dilation and strain ellipse (magnitude and orientation) of the upper surface digitalized with photogrammetry. Fold classification of Lisle & Toimil (2007): synclastic antiform (G>0, M>0), anticlastic antiform (G<0, M>0), synclastic synform (G>0, M<0) and anticlastic synform (G<0, M<0). Gaussian (G) and mean normal curvatures (M) are also displayed along with strike and dip maps; B) Dilation, strain ellipse and curvature patterns of the upper, middle and lower surfaces digitalized using cross section CT images. The strain ellipse gives us extra information about anisotropy. As observed before, upper and lower surfaces display opposite results. The same applies for the orientation 74 4 Surface restoration methods using paleomagnetism In this chapter, we introduce the two surface restoration methods that leverage paleomagnetic vectors as a primary reference. They are valid for complexly folded structures. The first one is a simple geometric approach based on the piecewise restoration of a triangulated surface (Gratier et al., 1991 and Gratier and Gullier, 1993) into which paleomagnetic variables can be easily incorporated (Ramón et al., 2012). The surface is modeled by a mesh and the method starts from a pin-element. Triangles are laid flat and then fitted together to minimize distances between common vertices and paleomagnetic error. However, this first approach, as it will show later, is sensitive to the meshing and particularly to the pin-element. The second one is based on a parametric approach, whereby a curvilinear coordinate system is computed on the folded surface by numerical optimization. We use paleomagnetic data to define constraints for the computation of this frame, which significantly increases the robustness of the restoration method. As in any restoration method, we need to establish some reasonable initial assumptions. We assume that: 1) Layers are horizontal in the undeformed surface; horizontality is the basic assumption in all restoration methods. 2) The folded surface is developable; it has been transformed by preserving angles, lengths and areas, so the Gaussian curvature is constant and zero. The method is also valid for globally developable surfaces, like those derived from rock volumes that have undergone flexural folding as described by Ramsay (1967), i.e., flexural flow on flanks or tangential longitudinal strain at hinges. In this case, total volume is assumed to be constant. 3) An even distribution of paleomagnetic vectors characterizes the folded surface. These vectors are primary (recorded at the time of deposition) and behave as passive markers during the deformation. In the analogs, these premises were clearly established, in the real databases, both the local and the paleomagnetic reference vectors have to be reliable in the sense used by Van der Voo (1990) and Pueyo (2010). 81 4.1 Piecewise approach The starting point of this approach is the UNFOLD method developed by Gratier et al. (1991) and by Gratier and Guillier (1993), as it is a simple geometric approach into which paleomagnetic variables can be easily incorporated (Ramón et al., 2012). We also considered an alternative method (Williams et al., 1997) that minimized internal deformation instead of distances, preserving the triangle’s area and allowing the modification of triangle’s shape. But the more proper method to introduce the paleomagnetic constrain is the one defined by Gratier et al. (1991) because triangles rotate rigidly. The UNFOLD method, like most flexural restoration tools, requires the first two of the previous assumptions plus the paleomagnetic restrictions, which can be easily integrated. The horizon is defined as a mesh of triangular elements that are first laid flat, and then rearranged (translated and rotated) to minimize distances between neighbors. The main change on adding paleomagnetic data is that the rotation is not free when minimizing the distances, as it is constrained by the paleomagnetic reference vector (Fig. 4.1). The software (Pmag3Drest) has been developed using Matlab1. Figure 4.1: The concept: paleomagnetism as an additional tool in restoration methods. The surface is rotated to the horizontal with the bedding plane and then it is rotated around its vertical axis to fit with its paleomagnetic reference vector. The method involves the following sequence of eight steps (Fig. 4.2): 1) Surface definition. A set of points with Cartesian coordinates describes the folded surface and determines the nodes of the rigid triangular elements of the mesh. To build the mesh, Delaunay (1934) or regular triangulation (Hjelle and Dæhlen, 2006) can be used. The most common triangulation is Delaunay’s one which states that the 1 www.mathworks.com 82 circumscribed circle of every triangle should not contain inside it any point of the mesh. It maximizes the minimum angle of all the angles of the triangles, which, as a practical result, avoids obtaining triangles with very small angles. However, in the analog models, we use the regular triangulation because it is much easier to obtain using commercial off-the-shelf wire meshes or by printing, in order to set the reference grid. In any case, we are going to evaluate the influence of the meshing type in the next chapter. 2) Incorporation of paleomagnetism. The method adds paleomagnetic vectors in triangular elements where paleomagnetic data is available, only considering the magnetic declination (the so-called horizontal component). The vector passes through the barycenter of the triangle (Fig. 4.3). The accuracy of the data, α95 (Fisher, 1953), can be related to the local paleomagnetic vector. In the case of sparse or poorly distributed vectors, it is possible to interpolate paleomagnetic data to all triangles of the surface using the algorithm described in Chapter 6. 3) Flattening. Each triangle is automatically laid flat to form a horizontal surface, by horizontal rotation about its strike axis. In the case of overturned beds, we treat the stratigraphic polarity as a vector in each element. * Pin-element definition. Following 2D restoration techniques, the unfolding restoration tools select a pin-element (also called the seed in the parametric approach), which is used as a fixed reference position in a restored model. Two types of pinelements can be defined for a surface: a pin-point is the first fixed triangle from where we start to unfold the surface, while the pin-line is the first row of triangles with fixed barycenters. The pin-element needs to be chosen with geological meaning because restoration is highly dependent on it, as it will be shown in next chapter. The pin is usually placed in the undeformed area of the surface which is considered not to have undergone deformation (i.e. the foreland). 4) Vertical-Axis Rotation. The paleomagnetic vector from the pin-element must converge with the paleomagnetic reference vector. Significant rotations are very unlikely since the pin-element should be chosen in the stable (undeformed) portion of the horizon. However, if there are substantial rotations, we apply an equal rigid-body rotation to the horizon at this stage. 83 5) Translation and rotation. Each triangle is translated and rotated in order to fit its neighbors using the method of least-squares. We minimize distances between shared vertices, bearing in mind the paleomagnetic reference. This step starts at the pinelement. If an element has paleomagnetic data, the related rotation is constrained by the paleomagnetic vector (Fig. 4.3), with a number of degrees of freedom that is determined by the α95 value. If the initial (undeformed) surface were completely developable (without any deformation), the restoration process would end at this step. 6) Iterating. The translation and rotation process is iterated a certain number of times, or for as long as the total distance error remains below a threshold, ∑∑ =MDe . In this expression, e is the error, D the sum of distances between the vertices of each triangle and the triangular hole defined by its neighbors and M the sum of the medians of each triangle (Fig. 4.3). This step is especially important in the original UNFOLD method, where paleomagnetism is not considered (Gratier et. al., 1991) and rotation is free. When rotation is constrained this step is unnecessary. 7) Welding. After the iterative translation and rotation, the surface becomes discontinuous, with gaps and overlaps, assuming that the deformed surface that has been restored was not completely developable. In order to obtain a continuous surface, this step involves joining or welding the shared vertices of neighboring triangles through an average value, allowing internal deformation to take place. 8) Optimization process. This step is only performed when the restoration uses paleomagnetic data. At this point of the restoration process, the maximum paleomagnetic error (the difference between the local and the reference paleomagnetic vectors) is less than the α95 angle, and we may want to sacrifice accuracy in favor of area preservation. The vertices of the triangles are randomly modified to minimize a potential function (Eq. [1]), following the simulated annealing method (Kirkpatrick et al., 1983; Press et al., 1992). The potential function (U) includes the paleomagnetic error (1cos[ref-pmagi]) and the internal deformation (dilation in terms of area variation: [area0-areai]/areai) with specific weights for each term (A and B). [] ∑∑ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛− +−−= ii i ii iarea areaarea BpmagrefAU 0 )cos(1 [1] This step was proposed for the restoration with paleomagnetism although practical experiments show that this step is avoidable because we lose primary information. The 84 optimization algorithm tends to distribute the error across the entire surface. It is useful to obtain a smoother surface, when we are able to ensure there is no deformation. On the other hand, it is better to omit this last step when we are trying to identify a possible deformation, as it causes the model to lose information; the deformation always becomes weaker after the optimization (Fig. 4.4). Although the optimization process has been considered, we end the restoration at step 7 in the simulations shown later. Figure 4.2: The method, step by step. 1) Surface definition. 2) Incorporation of paleomagnetic data. 3) Flattening. 4) Vertical axis rotation (unneeded in this example because paleomagnetic data of the pinelement fits with the reference). 5) Translation and rotation fitting. 6) Iterating process of step 5 with same results. 7) Welding of common vertices. 8) Optional optimization process. 85 Figure 4.3: Rotation of an individual triangle with and without using the paleomagnetic vector. c0 is the barycenter of the triangle. pmag is the paleomagnetic vector and ref its reference. Figure 4.4: Dilation of the restored surface before and after the optimization step. Due to the optimization, the dilation decreases but the expected pattern disappears. 4.1.1 Formulation We now detail the particular equations used in the code: For a triangle defined by its vertices vi, its barycenter is 3 321 0 vvv c++ = and its normal vector is 0201 cvcvN × = . The rotation matrix for the flattening step v is: 00 ')( ccvR ii +−⋅= 86 ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ +−⋅⋅+−⋅⋅⋅−−⋅⋅ ⋅−−⋅⋅+−⋅⋅+−⋅⋅ ⋅+−⋅⋅⋅−−⋅⋅+−⋅ = ψψψψψψ ψψψψψψ ψψψψψψ cos)cos1(sin)cos1(sin)cos1( sin)cos1(cos)cos1(sin)cos1( sin)cos1(sin)cos1(cos)cos1( 2 2 2 zxyzyxz xzyyzxy yzxzyxx R where xyz is the rotation axis and ψ the magnitude. As the rotation axis is horizontal it is possible to write the rotation matrix as: () () ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ΨΨ⋅Ψ⋅− Ψ⋅−Ψ+Ψ−⋅Ψ−⋅⋅ Ψ⋅Ψ−⋅⋅Ψ+Ψ−⋅ = coscos coscos)cos1()cos1(cos )cos1(coscos)cos1(cos 2 2 sensensen sensensen sensensen R φφ φφφφ φφφφ where )(90 90)( Ninc Ndec −= − = ψ φ If the surface is inverted: ψ ψ − = 180 inv The fitting process of translation and rotation is an iterative process: 1) Fix pin-element triangles. 2) Seek for neighbourhoods. 3) If pmag is defined Æ Rotate to fit with its reference: pmagref −= α , , ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛− = αα αα cossin sincos R00 ')( ccvRv ii +−⋅= 4) Translation and rotation. Objective: minimize distance between vertices. Distance: ()() ∑∑ == −+−== 3 1 22 3 1 ),( i ibiibi i ibi yyxxvvdistd Minimizing distance: 0,0,0 = ∂ ∂ = ∂ ∂ = ∂ ∂ α d y d x d New vertices: 00 ')( cossin sincos cntranslatiocvv ii ++−⋅ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛− = αα αα • 2 vertices fixed Translation when 2 vertices are fixed (rotation is determined by pmag): 2/)( 2211 abab vvvvntranslatio −+−= where a is the initial triangle and b the fixed. Translation and rotation with 2 vertices of the triangles fixed (only if pmag is not defined or if it has some error allowed; 0 95 > α ): 87 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⋅⋅⋅⋅⋅⋅⋅⋅ ⋅⋅⋅⋅⋅⋅⋅⋅ = b1a2b2a1b1a2b2a1b2a2b2a2b1a1b1a1 b1a2b2a1b1a2b2a1b2a2b2a2b2a1b1a1 yy-yy-xx-xx-yy+xx+yy+xx yx-yx-xy+xy+yx+xy-yx+xytana α If 9595 α α α α =→> , ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⋅⋅⋅ ⋅⋅⋅ =)cos)y+(y-sin)x+(x-y+(y1/2 )sin)y+(y+cos)x+(x-x+(x1/2 a2a1a2a1b2b1 a2a1a2a1b2b1 αα αα ntranslatio • 3 vertices fixed Translation with 3 vertices of the triangle fixed. The triangle must have the same centroid than the hole made by the mean vertices of its neighbours: 00 ab ccntranslatio −= Only rotation with 3 vertices of the triangles fixed: ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⋅+⋅+⋅+⋅+⋅+⋅ ⋅−⋅+⋅−⋅+⋅−⋅ = 333322221111 333322221111 tan babababababa babababababa yyxxyyxxyyxx yxxyyxxyyxxy a α 5) After fitting each triangle, come back to the point 2. The optional optimization step is based on the Simulated Annealing method. Nodes are moved randomly to decrease a potential function U (Eq. [1]). The change is accepted if and if 01 UU < () ae Tk UU > ⋅ − −01 (some local worsen is allowed in order to find the global minimum and not only a local one). The temperature (T) decreases with the number of iterations, at the end only changes that really decreases U are allowed α ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + −= 1 1 0M m TT where m is the actual iteration and M the total number of iterations. By default: 1,1,02.0,9.0 0 = === α Tka 4.2 Parametric approach This proposed alternative approach incorporates paleomagnetic information in an unfolding algorithm based on the parameterization of the surface, a method also valid for surfaces folded under flexural conditions (Ramón et al., in review). We have incorporated this constraint in the gOcad (by Paradigm2) code (Massot, 2002). 2 http://www.pdgm.com 88 Any surface S ( ) can be projected onto a map C ( ) using a parametric representation (Fig. 4.4). This so-called surface parameterization problem has been extensively studied in the computer graphics community. We summarize the main ideas below and refer the reader to Botsch et al. (2010) and Floater and Hormann (2005) for further details and discussions, and Mallet (2002), Moretti et al. (2005) and Moretti (2008) for more details about the application of the free boundary parameterization to surface restoration. Each point in C defined by the coordinates u=(u,v) has a unique associated image in S x(u) with coordinates (x,y,z): 3 ℜ2 ℜ S z y x C v u∈ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎯⎯⎯← ⎯⎯→⎯ ∈ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =xu xu ux )( )( [2] For a single surface there are many possible parametric representations. However, for any parameterization it is possible to define the metric tensor G(u,v) which defines the intrinsic properties of the surface. The main assumption of most restoration algorithms is the horizontality of the initial surface. Thus, the restoration problem is equivalent to look for a particular parameterization u(x). Moreover, we focus on surfaces folded under flexural conditions that lead to globally developable surfaces (Gaussian curvature equal to zero almost everywhere). This assumption implies the preservation of lengths and angles in the folded surface. The paleo-geographic coordinate system (u,v) can be freely and arbitrarily chosen and because of simplicity it is assumed to be rectilinear and orthonormal. Because of the principle of minimum deformation, the metric tensor G associated with the local parameterization must be close to the unit tensor (Eq. [3]). () ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =≅→≅−⋅= 10 01 ),(0),(2/1 IvuGIvuG ε [3] From a given 3D surface, the restoration problem can be addressed by computing two coordinates (u,v) on the surface so that: 1 0 =∇=∇ =∇⋅∇ vu vu [4] 89 The first constraint means that both gradients are orthogonal and the second that they have the same norm (Fig. 4.4). These two conformality constraints can be discretized on the surface in order to compute the parameterization, e.g. see Lévy and Mallet, (1999) or Levy et al. (2002). In this work, we started with a parameterization method developed by Massot (2002), which starts from a pin-element on the surface where two directions and are chosen so as to be mutually orthogonal and orthogonal to the local normal vector n0 (Fig. 4.4B). This method starts by propagating this local frame to all the nodes of the surface. For this, the frame is rotated along an axis l given by the cross product of the normal at the pin-point n0 and the normal at the current point n1. (Figs. 4.4B&C). Last, the u and v coordinates are integrated numerically from the pin-point coordinate using a least-squares method so as to honor the input gradients at all locations (Fig. 4.4D), see details in Massot, (2002) and Mallet (2002). 0 u∇0 v∇ Although this method often provides good results in moderately deformed areas, it remains quite sensitive to the initial solution and fails for complexly folded structures. We improve this restoration approach by using paleomagnetic data. The idea is to use the paleomagnetic values defined in all points of the surface as the initial gradient values . The gradient is a vector contained in the surface, thus, the initial paleomagnetic vector is rotated u∇ to have null inclination and so to be embedded in the surface. Paleomagnetic data is known in specific sites and then interpolated to the entire surface as described in Chapter 6. 90 Figure 5.2: Curvature analysis and real deformation of the San Marzal model. Properties are displayed in the initial horizontal surface. Dashed black lines divide the maximum dilation areas. A) Mean normal curvature (M=(k1+k2)/2). B) Gaussian curvature (G=k1*k2). C) Fold classification described by Lisle and Toimil (2007): 1.synclastic antiform (M>0, G>0), 2.anticlastic antiform (M>0, G<0), 3.synclastic synform (M<0, G>0), 4.anticlastic synform (M<0, G<0). D) Real dilation (dreal=(Afolded-Ainitial)/Ainitial). E) Anisotropy of the strain ellipse; strain ratio: major axis / minor axis. F) Strain orientation. Major axis orientation. Axis displayed in those triangles with meaningful strain ratio (>1.07). 97 Figure 5.3: Curved fold. A) Analog model based on the Balzes Anticline. B) Digitalized model with values of normal curvature (M) displayed. The small picture shows the inverted area. Figure 5.4: Curvature analysis and real deformation of the Balzes model. Dashed black lines divide the maximum dilation areas and dashed blue lines the maximum strain ratio areas. A) Mean normal curvature. B) Gaussian curvature. C) Fold classification: 1.synclastic antiform, 2.anticlastic antiform, 3.synclastic synform, 4.anticlastic synform. D) Real dilation. E) Strain ration. F) Strain orientation. Major axis orientation. Axis displayed in those triangles with meaningful strain ratio (>1.07). 98 A similar analysis is done for the second analog model; the Balzes Anticline (Fig. 5.4). Again, maximum dilation areas correspond with those of maximum curvature, all located at the hinges of the folds where tangential-longitudinal strain has taken place. However, in this case, the maximum strain ratio has a bit different pattern and it is possible to associate it with the Gaussian curvature. The most deformed areas are those with maximum Gaussian curvature (no matter the sign). The orientation of the major axis of the strain ellipse seems to be perpendicular to the fold axis orientation in the outer hinges (as expectable from the Gairola’s [1967] models) where there is dilation and parallel to the fold axis in the inner hinges where there is contraction. 5.2 Restoration with and without paleomagnetism After this round through the real deformation parameters we are going to test the incorporation of the paleomagnetic constraint into the restoration methods. We are going to unfold the surfaces with the piecewise approach using and not using paleomagnetism (Pmag3DRest and UNFOLD methods respectively) as well as with the parametric approach (gOcad method with and without the paleomagnetic constraint). The piecewise restoration is quite sensitive to the pin-element selection, as we will see in Section 5.4, but for our current purpose we stick to the traditional approach ofplacing a pin-element in the undeformed area, although not in the border. The orientation of the paleomagnetic vector also has some influence in the results; for this example we try to select one oblique to the main structures that better registers the deformation, but of course this depends on the real geology setting in natural cases. Eventually, in these simulations we consider paleomagnetism defined in all triangles. We first analyze the piecewise restoration of the San Marzal model using and not using paleomagnetism (Fig. 5.5, Table 5.1). Restoration without paleomagnetism seems to be pretty good in terms of mean dilation (dmean without pmag < dmean with pmag) but restoration with paleomagnetism better locates real deformation (e with pmag < e without pmag). The restoration without paleomagnetism achieves the restored surface with minimum deformation because triangles are free to find the best fitting. However, the paleomagnetic constraint yields a restored surface more similar to the initial one. Dilation values of the restored surface with paleomagnetism are mostly coherent with 99 the real dilation in the expected areas: expansion in the anticline and compression in the synclines. On the other hand, absolute values are always smaller than those expected (dmean < dreal_mean, e >> 0) because of the restoration method itself which try to minimize deformation. Additionally, the anisotropy of deformation (strain ellipse) is coherent for the restoration with paleomagnetism and has less meaning for the restoration without paleomagnetism. The compression area with maximum strain ratio coincides with the real one and shares the same orientation. Moreover, the shape of the surface restored with paleomagnetism is closer to the expected rectangle. This is a good example where piecewise restoration works well and paleomagnetism provides additional information about the real deformation of the surface. Let us compare now this results with the second method based on the parameterization. The parametric approach without the incorporation of paleomagnetism is unable to reach a proper restoration of the surface (Fig. 5.6B) and mean dilation values are to high (dmean ≈ 10% [Table 5.1]). We can say that this restoration method is only valid for smoothly folded structures and not for complex folds such as this example. This method is too dependant on the initial solution and for this reason the usage of paleomagnetism is crucial to define it (particularly in non-cylindrical and non-coaxial structures). The shape of the restored surface with paleomagnetism almost perfectly matches the expected one and orthogonality is preserved in the restored state. The distribution of maximum dilation values roughly corresponds with those expected, although the signs do not match in this case (restored surface does not show the contraction area). On the other hand, maximum strain orientation approximately coincides with the expected one. The restoration method assumes isometric parameterization based on developable surfaces folded under flexural conditions but we have used it to restore globally developable surfaces with small deformation in complex folded areas. Therefore, the unfolding algorithm minimizes the deformation showing dilation patterns similar to those real, but it could never predict the real deformation. Because of that, we recommend the usage of several restoration methods in order to compare results. In this case for example, the strain orientation of the compression area with a maximum strain ratio is trustable because almost coincides with both techniques. 100 Figure 5.5: Piecewise restoration of the San Marzal model (regular triangulation, original mesh density). Dashed black lines are useful to compare the restored surface with the initial one and real values of deformation (Fig. 5.2). The pin-element is highlighted with a red circle. A) Restoration using paleomagnetism (Pmag3DRest method). The red arrow is the paleomagnetic orientation. B) Restoration without paleomagnetism (UNFOLD method [Gratier et al., 1991]). 1) Dilation (d=(AfoldedArestored)/Arestored). 2) Strain ratio: major axis / minor axis. 3) Strain orientation: major axis displayed in those triangles with meaningful strain ratio (>1.07). 101 Figure 5.6: Parametric restoration of the San Marzal model. Dashed black lines are useful to compare the restored surface with the initial one and real values of deformation (Fig. 5.2). A) Restoration using paleomagnetism: 1) dilation, 2) strain ratio, 3) strain orientation. B) Restoration without paleomagnetism with dilation values. SAN MARZAL (dreal_mean=0.0433) Piecewise restoration Parametric restoration with pmag without with pmag without e 0.0392 0.0442 0.0459 0.0964 dmean 0.0174 0.0086 0.0164 0.0925 Table 5.1: Mean dilation values for the restoration of San Marzal. Mean error between real and retrodilation in absolute value: ∑−= ddNe real 1 (Eq. [8]) Mean dilation in absolute value: ∑ =dNdmean 1 (Eq. [6]). Mean real dilation in absolute value: dreal_mean. The next example (curved fold based on the Balzes Anticline) is more complicated than the previous because it has more deformation (higher dreal_mean) and the surface is less developable (higher Gaussian curvature). Moreover, the meshing is denser, hampering the piecewise restoration. 102 Piecewise restoration (Fig. 5.7, Table 5.2) is not as accurate for this example and the benefits of paleomagnetism are not as clear (e & dmean without pmag < e & dmean with pmag). Dilation increases towards the opposite part of the pin-element due to the propagation of errors. Moreover, the dilation values seem to be more related to the areas of higher Gaussian curvature than with the real dilation, particularly in the interference syncline between both anticlines. Figure 5.7: Piecewise restoration of the Balzes model. Dashed lines are useful to compare the restored surface with the initial one and real values of deformation (Fig. 5.4). The pin-element is highlighted with a red circle. A) Restoration using paleomagnetism (Pmag3DRest). The red arrow is the paleomagnetic orientation. B) Restoration without paleomagnetism (UNFOLD). 1) Dilation. 2) Strain ratio. 3) Strain orientation. 103 Figure 5.8: Parametric restoration of the Balzes model. Dashed lines are useful to compare the restored surface with the initial one and real values of deformation (Fig. 5.4). A) Restoration using paleomagnetism: 1) dilation, 2) strain ratio, 3) strain orientation. B) Restoration without paleomagnetism with dilation values. In this situation, it seems more convenient to unfold the surface with the second restoration method, the parametric approach, because the results are really encouraging (Fig. 5.8). As in the example of San Marzal, the method does not work well without the paleomagnetic constraint. Thus, we focus on the restoration with paleomagnetism. Dilation and strain ratio patterns are fairly similar to those expected. Strain orientation is more confuse this time. BALZES (dreal_mean=0.0989) Piecewise restoration Parametric restoration with pmag without with pmag without e 0.1087 0.0988 0.0856 0.1166 dmean 0.0415 0.0240 0.0228 0.0572 Table 5.2: Mean dilation values for the restoration of Balzes. 104 5.3 Mesh sensitivity In discrete computer-based models some algorithms may lead to artifacts when surface sampling changes. Therefore, we check the robustness of the method towards variations in density and structure of the mesh (Ramón et al., in review). 5.3.1 Mesh density sensitivity The impact of mesh density is evaluated by decimating the initial point set defined by the nodes of the reference grid and preserving the meshing method and boundary of the surface. We decimate the initial nodes by a factor of two for rows and columns. The San Marzal model initially has 722 triangles and 400 nodes and after decimation: 200 triangles and 121 nodes. The Balzes model has 3570 triangles and 1872 nodes and after decimation: 936 triangles and 513 nodes. The impact of mesh density for the simple geometry model (San Marzal) is only evident when analyzing the values of dilation for the restored surface because it has a proper rectangular shape (Fig. 5.9). With a higher mesh density, dilation in the restored surface is more localized along the hinges of the fold, whereas the lower mesh density model has a greater dispersion of dilation, partly due to the difference in triangle size (Fig. 5.9 vs. Fig. 5.5 & Fig. 5.6). However, differences in mean dilation between low and high density models are not great (Table 5.3 vs. Table 5.1), indicating that mesh density has a limited impact on restoration results. 105 Figure 5.9: Mesh density analysis of San Marzal (regular triangulation). A) Decimated initial surface (with boundary preservation). B) Restored surface using the piecewise approach with paleomagnetism. Dilation values displayed. C) Piecewise restoration without paleomagnetism. D) Parametric restoration with paleomagnetism. E) Parametric restoration without paleomagnetism. Similar considerations can be done for the second model of a curved fold (Balzes) (Fig. 5.10). Areas of maximum and minimum dilation obtained with both restoration methods are equivalent for the surfaces with more and less mesh density (Fig. 5.10 vs. Fig. 5.7 & Fig. 5.8). However, in this case there is an additional error produced by the influence of the border, particularly in the neighborhood of the pin-element where the size of triangles is smaller in order to preserve the boundary of the initial surface. This error causes an increase on the mean dilation value (Table 5.3 vs. Table 5.2). 106 Figure 5.13 (continued): Pin-element analysis, piecewise restoration starting from different points. 1) San Marzal model. 2) Balzes model. A) Restoration with paleomagnetism. B) Restoration without paleomagnetism. Dilation values displayed. The restoration of the Balzes model is more uncertain in any case. Dilation patterns of restored surfaces highly vary as the pin-element placement changes. This difference among solutions raises the question of whether the deformation obtained with the restoration represents the real one or not. In any case, the restoration with 113 paleomagnetism clearly improves the solution yielding the proper rectangular shape of the initial surface. Another consideration we can make by observing the results (particularly from the Balzes model) is that generally, the deformation concentrates far from the pin-element. Therefore, if we want to use the restoration to predict real deformation, we should avoid placing the pin-element in a deformed area (as in the first simulation of Fig. 5.13). Figure 5.14: Pin-element analysis for the parametric restoration without paleomagnetism. 1) San Marzal model. 2) Balzes model. Dilation values displayed. The parametric restoration is sensitive to this variable only when restoring without paleomagnetism, on the other hand, only useful to restore simpler structures. The main difference is the selection of the pin-element inside or outside the main structure (Fig. 114 5.14): a pin-element located inside the main structure where deformation has taken place seems to be a better choice in order to get a restored surface more similar to the initial one. 5.5 Multi-surface restoration The restoration methods developed in this work unfolds single surfaces and do not consider volumes. However, multi-surface restoration is the first step for a multi-map volume restoration and it may be a good approximation to the real 3D problem. Figure 5.15: Multi-surface restoration of the Balzes model. A) Initial surfaces (upper and lower) with real dilation values. B) Restored surfaces with retro-dilation values. Note: the range of dilation is not the same for initial and restored surfaces, being higher for the first. Thanks to the CT scanner we are able to reconstruct internal surfaces of the models and we can characterize completely the Balzes, although the resolution is worse than in the case of photogrammetric reconstruction. The restoration method we are going to apply to the model obtained with the CT scanner is the parametric approach with orientation of paleomagnetic data of pmag3 because it is the one that has shown better results. We restore separately the upper and lower surfaces of the model (just the most 115 representatives) and compare them with the expected results. Real dilation values (already analyzed in Section 3.2.3 Fig.3.17) show extension in the outer parts of the anticline and compression in the inner parts, in accordance with the tangentiallongitudinal deformation described by Ramsay (1977). Restored surfaces (Fig. 5.15) display opposite signs of dilation which suggest the real behavior of the volume: extension in the outer parts of the anticline and compression in the inner parts. Dilation is higher in the real case than in the restoration, but the consistency in the difference between both surfaces (opposite dilation signs) is quite encouraging. 5.5.1 Comparison with a 3D restoration method We may think that multi-surface restoration is far from recent real 3D restoration methods based on geomechanical approaches. However, for highly folded surfaces as those of our examples, these methods assume elasticity of materials and present several limitations we analyze in this Section We used Dynel3D (IGEOSS1; Maerten & Maerten, 2006) to restore the Balzes model. The stratigraphic units in Dynel3D are discretized with tetrahedral elements that are assigned elastic properties. The tetrahedral elements deform elastically in response to constraints such as applied and/or internal forces, displacements and interface contacts. Dynel3D uses an iterative, explicit solver that preserves mass and allows forces to be transmitted from node to node through the entire system until equilibrium is reached. One of the characteristics of the finite element methods (FEM) as this one, it is that requires strong boundary conditions, treating the geological body as a separate individual system delimited by a volume within which the calculations take place. A closed boundary or bounding box must be defined prior to the volume generation (Fig. 5.16A). On the other hand, the physical-based restoration algorithm run in Dynel3D needs several rock mechanical properties to be set up. As these mechanical properties vary with lithology, we set the typology of a sandstone (Young's modulus: 2.2E10Pa, Poisson's ratio: 0.24, density: 2480kg/m3). The method also permits to model the 1 http://www.igeoss.com 116 behavior of the contact between units. Since the two layers of the model represent the same stratigraphic unit, we have blocked them in Dynel3D. Figure 5.16: Volumetric restoration of the Balzes model with Dynel3D. A) Initial model with its boundary box. B) Restored upper layer. C) Restored lower layer. Volumetric strain displayed: E1+E2+E3. With these initial settings, we unfold the layers assuming horizontality of initial state. The properties calculated are different than the used in our methods, but we can qualitatively compare results. We display the volumetric strain which is the sum of all main strains. The restoration method is unable to undo the deformation process because it can reach the initial rectangular shape of the surface (Figs. 5.16B&C). Strain consistently accumulates in the hinge anticlines of the upper layer and hinge synclines of the lower layer. However this strain is always positive, indicating than the volume needs to expand to reach the folded surface, and therefore, an unreal change of volume has occurred. 5.6 Conclusions In this chapter, we have evaluated the goodness of the piecewise and parametric restoration methods using paleomagnetism, which is defined in all points of the surface. The paleomagnetic constraint improves the resultant surface in most cases achieving in general more similar deformation patterns to those expected. Particularly, we can highlight the importance of paleomagnetism when the methods present higher limitations: 1) the strong sensitivity to the pin-element of the piecewise restoration; and even more important, 2) the impossibility of the parametric restoration to achieve the expected solution for strongly folded surfaces (Table 5.5). 117 Additionally, we recommend the multi-surface restoration as previous step to real 3D restoration. If horizons are correctly modeled and the surface restoration method applied works properly for one surface; the multi-surface restoration is able to indicate differences between surfaces offering volumetric information. We also emphasize the importance of CT analogs to test 3D restoration methods. Piecewise rest. Parametric rest. Restoration Rest. with pmag better locates deformation patterns Rest. w/o pmag is invalid for complex structures Sensitivity to meshing Dilation patterns are equivalent Dilation patterns are equivalent Sensitivity to pin-element Too sensitive (although pmag improves the results) Insensitive with pmag Multi-surface restoration Good 3D approximation Table 5.5: Summarizing table. 118 6 Checking the consistency of the restoration methods: paleomagnetic sensitivity In this chapter we analyze the influence of the initial paleomagnetic dataset on the results of the restoration methods. Specifically we analyze: - The influence of the initial orientation of paleomagnetic vectors. - The sensitivity to the error degree inherent to the definition of paleomagnetism for a given site. - The sensitivity to data availability and distribution, since in a real scenario paleomagnetic data is not defined in all points of the surface, being only available in scattered sites. Additionally, we propose an interpolation algorithm to extend sparse data to the whole surface. Then, we analyze the usage of scattered and interpolated data. 6.1 Sensitivity to the initial paleomagnetic orientation In this section we analyze the effect on the orientation of the paleomagnetic vectors. In a real case we can not modify this parameter but we need to know if the orientation of initial paleomagnetic data conditions the result. Paleomagnetism in the analog models is defined by the sides of the triangles, and therefore we have three possible initial datasets. We restore the two analog models with the two restoration methods using the three paleomagnetic datasets (Fig. 6.1). The San Marzal model presents two preferred orientations with lower mean dilation and error values (Table 6.1) and better location of maximum dilation areas: pmag1 (dataset used in previous examples) and pmag3 although the observed dilation is always smaller than the expected one. These orientations are oblique to the main structure, while pmag2 is rather parallel to the fold axis. Although with this single observation is not feasible to generalize the result. The restoration of the Balzes model with the piecewise approach is more uncertain than the restoration with the parametric approach and, accordingly, we especially focus on the second method. The three restored surfaces present positive dilation in the main 119 anticline and negative dilation in the main syncline between both folds. Main differences appear in areas of maximum Gaussian curvature and at the boundary. Between the three simulations, the one restored with pmag3 has the minimum mean dilation error (e, Table 6.1) and presents the dilation pattern closest to the expected. This is the dataset used in previous examples, and again, is the most oblique to both folds. These results show that the methods examined are quite robust to the paleomagnetic data orientation. Although there are some variations in the dilation patterns, they are much less significant than the ones we observed for other factors, like for example the pin-element variation. That is quite reassuring, because, as mentioned before, in a real case we can not vary the paleomagnetic data. 120 Figure 6.1 (continued): Paleomagnetic orientation analysis. Paleomagnetism is defined with the different sides of the triangles: pmag1, pmag2, pmag3. 1) San Marzal model. Previous examples used pmag1. 2) Balzes model. Previous examples used pmag3. A) Piecewise restoration with paleomagnetism. B) Parametric restoration with paleomagnetism. Dilation values displayed. pmag1 pmag2 pmag3 e 0.0392 0.0519 0.0391 Piecewise rest. dmean 0.0174 0.0191 0.0123 e 0.0459 0.0473 0.0451 San Marzal model Parametric rest. dmean 0.0164 0.0186 0.0169 e 0.0993 0.1051 0.1087 Piecewise rest. dmean 0.0460 0.0415 0.0415 e 0.0924 0.0952 0.0856 Balzes model Parametric rest. dmean 0.0198 0.0177 0.0228 Table 6.1: Mean error and dilation values for restored surfaces using different initial datasets. 121 6.2 Sensitivity to the paleomagnetic accuracy Paleomagnetic data must be a reliable record of the magnetic field at the time of rock formation (Van der Voo, 1990; Pueyo, 2010) and it is always defined with a given accuracy (α95). In this section we analyze the effects of the paleomagnetic error degree. We define two initial datasets: 1) with a random error ≤5º with a confidence level of 0.95 (α95=15º) and 2) with a random error ≤15º with a confidence level of 0.95 (α95=15º). We restore the two analog models with the two restoration methods using both paleomagnetic datasets (Fig. 6.2). For the piecewise restoration there are two options: 1) use the paleomagnetic data as hard data although it is not completely accurate (Fig. 6.2B) and 2) use it with its angular variation (α95) allowing a free rotation of triangles according to that (Fig. 6.2C). This second option needs the iterative step used in the restoration without paleomagnetism (step 6 described in Section 4.1). In order to quantify the difference between the restored surface using and not using accurate data, we measure the mean error dilation between the surface restored from accurate paleomagnetic data and the surface restored (with the same method) from inaccurate paleomagnetic data (erest=mean(|drest0-drestAlfa95|, Table 6.2). In this case we do not compare the real dilation with the restored one (e) as we did in previous cases, but the dilation of the restored surface using ideal and real initial datasets (erest). For the initial dataset defined with an error minor than 5º and for both models, the piecewise restoration that uses initial data as hard data, without allowing free rotation of triangles, (option 1) produces lower mean error dilation (erest) than the piecewise restoration that allows free rotation (option 2). This agrees with dilation patterns. It can be observed that the first option leads to similar restored surface patterns than those obtained with the piecewise restoration using paleomagnetic data (Fig. 6.2A vs. Figs. 5.5A & 5.7A) while the second option is similar to the restoration without paleomagnetism (Fig. 6.2B vs. Figs. 5.5B & 5.7B). In the same way, for the initial dataset with α95=15º, the dilation patterns of piecewise restored surface without free rotation (option 1 vs. option 2) are similar to the dilation patterns of the restoration with paleomagnetism and α95=0º, although the former have much more deformation. Because of that, the mean error dilation increases this time. Due to this deterioration of the results, we encourage to do an effort to achieve a reliable paleomagnetic dataset. We particularly propose in the appendices two ways to 122 Figure 6.4: Paleomagnetism interpolation. A) Rotation of the paleomagnetic vector to propagate it from triangle to triangle. B) Propagation order starting from two different initial sites. The next step computes the average of both paleomagnetic propagated datasets. C) The interpolated paleomagnetism is calculated using all points of the surface as initial data. We plot the difference in degrees between real and interpolated paleomagnetic data. D) Initial paleomagnetic sites location simulating a real case with sparse data. E) Difference between real and interpolated paleomagnetic data from the initial sites plotted in figure D. 129 6.5 Real-world case: sensitivity to accuracy and resolution Finally, we simulate a real scenario with inaccurate and scattered initial data that are interpolated. We reproduce the scattered datasets described in Section 6.3; 14 sites for the San Marzal model and 15 sites for the Balzes model (Fig. 6.3B). For these sites, we define two different initial paleomagnetic datasets, one with total accuracy (α95=0º) and the other with a random error ≤15º with a confidence level of 0.95 (α95=15º). These datasets are interpolated with the algorithm described in previous section. First of all, we need to find the error introduced by the interpolation algorithm. To do so, we calculate the difference between the real dataset with paleomagnetism defined in all points and the interpolated datasets from the scattered sites (both defined with and without total accuracy) (Fig. 6.5A). The simpler model has mean errors of 1.47º and 2.69º for the interpolated datasets of α95=0º and α95=15º respectively, while the more complicated model has mean errors of 2.92º and 3.54º. Again, error is never too big in comparison with the usual paleomagnetism resolution (α95 ≈ 10º). Secondly, we restore the surfaces using the interpolated datasets and compare the results with those obtained with paleomagnetism defined in all points (Figs. 5.5A, 5.6A, 5.7A and 5.8A). For the piecewise restoration we must clarify the following issue; in the translation and rotation step, we do not allow a free angular variation related to the α95 value, but use paleomagnetism as hard data. We do it this way because if we set a free variation of 15º (as α95) the restored surface is equivalent to the one restored without paleomagnetism as observed in Section 6.2. Observing the restoration of the San Marzal model we can say that the restoration method based on the parameterization of the surface becomes more robust. Results define with more accuracy the main deformation area in the anticline hinge (maximum dilation). Even more important is that dilation patterns are similar in all three situations. This observation is stronger in the Balzes model; although dilation patterns become diffuse and are closer to a Gaussian curvature than to the real dilation, they are more coherent for the parametric restoration than for the piecewise restoration. In order to quantify the difference between the restored surface using interpolated dataset and the initial one, we measure the mean error dilation between the surface 130 restored with paleomagnetism defined in all points and the surface restored with interpolated data (eRest=mean(|drest0-drestInterp|)). The method based on the parameterization has lower values (Table 6.4) in both examples. That means, in general terms, that parametric restoration is more stable than piecewise restoration. On the other hand, restored surfaces of the Balzes model are almost independent on the initial paleomagnetic dataset (α95=0º or α95=15º). The reason is that the paleomagnetism is averaged with the interpolated algorithm, and the interpolation of paleomagnetic values with more or less accuracy leads to similar results. This makes the method appropriate to be used in real cases. 131 Figure 6.5 (continued): Paleomagnetic accuracy and resolution analysis. Initial scattered paleomagnetic datasets already defined in Figure 6.3B with total accuracy (α95=0º) and with random error minor than 15º (α95=15º). 1) San Marzal model. 2) Balzes model. A) Paleomagnetic error (in degrees) between the real dataset with paleomagnetism defined in all points and the interpolated dataset from the scattered sites. B) Piecewise and C) parametric restoration. eRest Interp. α95=0º eRest Interp. α95=15º Piecewise rest. 0.0191 0.0272 San Marzal model Parametric rest. 0.0079 0.0113 Piecewise rest. 0.0681 0.0825 Balzes model Parametric rest. 0.0203 0.0208 Table 6.4: Mean error dilation between restored surface with paleomagnetism defined in all points and restored surface with interpolated data (erestInterp=mean(|drest0-drestInterp|)). 132 6.6 Conclusions Summarizing, the proposed restoration methods remain valid for a real case with scattered paleomagnetic sites, thanks to the interpolation algorithm. However, the parametric approach is more robust than the piecewise approach. On the other hand, we should not forget that the dilation patterns obtained with the restoration are just an aid to determine the real dilation, but they never provide the real dilation because they are conditioned for the restoration procedure. 133 134 7 Map-view restoration In this chapter we explore the usage of paleomagnetism in map-view restoration techniques. Map-view restorations based on balanced cross-sections, also known as palinspastic restorations, have been long used to understand 3D patterns in fold and thrust belts (Kay, 1945 and 1954). Together with structural, isolith and isopach maps, they are very useful in underground oil and gas exploration. Palinspastic maps display folded or faulted strata restored to their paleogeographical location before deformation took place (folding or faulting). Classic (see for example Dalhstrom, 1969), or even recent approaches (Price and Sears, 2000), are based on balanced cross sections (2D). To ensure the map-view reconstruction effectively considers the entire volume of the materials represented and, therefore, respects the lengths of lines and the thicknesses of individual layers. The term palinspastic is derived from the Greek palin meaning again, and spastikos meaning pulling (Allaby and Allaby, 1999). We could consider this approach as an early 3D restoration technique. The map-view restoration approach proposed in this chapter does not concern surface restoration of individual structures; here we propose a cartographic map view restoration applicable at local and regional scales. The main goal of this map-view restoration method is to unravel the vertical axis rotations. We first need to establish some assumptions (and derived limitations) as well as the aims and scope of the method: 1) Vertical Axis Rotations (VARs) may produce severe room problems in fold and thrust belts. Therefore, the removal of VARs at the regional scale is equivalent to the removal of the lateral gradient of shortening (in the sense used by Pueyo et al., 2004 and Sussman et al., 2012). This means that only parallel and constant translation will remain in the thrust front after the restoration of the map before the rotational period. The restored map will not represent a real image. Additional work should be done in the future to implement the combined restoration of VARs and translations due to cylindrical folding and thrusting. 2) The second assumption deals with the rotational time. Several evidences from the Pyrenees point to a narrow rotational window during the Eocene-Oligocene: LutetianPriabonian in the Southwestern (Pueyo, 2000; Pueyo et al., 2002; Mochales et al., 2012; 135 Rodríguez et al., 2013c; Muñoz et al., 2013) and in the Southeastern Pyrenees (Sussman et al., 2004), and similar ages both in the Northwestern (Mouleon basin) and in the Northeastern (Corbiers) (see works by Oliva et al. (2010) and Rouvier et al. (2012) respectively). Therefore, the map-view restoration will offer an imprecise snapshot of that moment, just before the out-of-plane movements began. 3) We also assume the geological map as a flat and horizontal surface. The final goal is to observe the structural trends before the rotational period. Within these three assumptions, our cartographic restoration technique aims to localize the expected gaps and overlaps between in relation to their cartographic location, and does not pretend to be a trustworthy reconstruction of the past. An additional advantage of our map-view reconstruction is its ability to identify the anisotropy related to the gaps and overlaps. Extension of surface restoration (2.5D) techniques by Audibert (1991) and by Rouby et al. (1993) using paleomagnetic vectors were developed by Arriagada et al., (2008) and applied to map-view (2D) in the Bolivian orocline. This technique, based on the least-square minimization of fault-boundaries blocks, assumes that the initial surface (geologic map) can be entirely divided by discontinuities in a set of discrete domains (heave maps). Although some interesting results can be derived from this technique, the problem is that many heave faults lack for geological expression or meaning, or at least they are scale-dependent in the sense used by Rouby et al., 2000, and this adds uncertainty to the restoration results. Our idea is somehow similar to the one developed by Arriagada et al. (2008), but in our approach, we divide the map in random triangles and not in blocks bounded by faults and heaves. In other words, we keep the continuity of the map. As we said, the map-view restoration only considers the coordinates xy, and structures are roughly assumed flat. We present two techniques based on the piecewise and parametric approaches. The first technique is equivalent to the fitting step of the piecewise restoration detailed in Section 4.1, which comes after the flattening step. Gaps and overlaps of the restored map (before the welding step) indicate areas were maximum deformation has taken place. After welding we can visualize the cartographic image in the restored map to approximate the initial location of the structures. This restoration technique is 136 essentially based on the rotation determined by paleomagnetic data: triangles with known data rotate according to this constraint, while the others fit minimizing distances between vertices. The second technique is the parametric restoration described in Section 4.2 in which the rotation data determines the gradient of one of the parameters. Whereas the drawback in the first technique is the pin-element location, the drawback in the second technique is that rotation data needs to be established everywhere. Again, paleomagnetism is the key-stone that allows quantifying vertical-axisrotations (VARs). However, VAR magnitudes can be now derived not only from primary vectors (like in the surface restoration) but also from secondary components (synfolding or postfolding) since they are still valid to quantify and date VARs (see Pyrenan examples in Oliva and Pueyo, 2007b). On the other hand, VARs are known in specific sites while the structural trend (or strike direction) can be easily determined in more points using geologic maps or field data. The proposed idea is to find some law to relate these two variables in order to have a large rotation dataset relevant for the restoration. To cope with this aim, we develop two distinct procedures: firstly, the use of specific vectors as the initial dataset (raw data), and secondly, leveraging strike vs. rotation relations, interpolated following the strike-VAR law. We apply this palinspastic restoration method at two different scales in two case studies, the Balzes Anticline (regional scale) and the South Central Pyrenees (tectonic scale). A vast paleomagnetic dataset is available in both cases (Rodríguez-Pintó et al., 2013c and López et al., 2008 respectively) 7.1 The method This palinspastic or map view restoration method starts from the cartographic map of the area of interest and the georeferenced rotation dataset, obtained from punctual paleomagnetic data [Section 7.1.1] or from paleomagnetic and structural data [Section 7.1.2]. The procedure is based on the piecewise approach as stated above, and therefore it involves a sequence of five steps: 1) Map triangulation. As in similar methods, the first step consists in the discretization of the area we want to restore. The cartographic map is meshed with a 137 Delaunay triangulation, targeting the homogeneity of the triangles. As discussed in the previous chapter (Section 5.3) the type of mesh hardly conditions the result, and therefore other approaches like regular triangulation can be employed. However, mesh density must be carefully chosen. It must be dense enough to represent all the structures and rotation data, and wide enough to minimize propagation errors. 2) Incorporation of rotation data and pin-element. The rotation vectors are added in the barycenter of triangles wherever they are known. Rotation data can be obtained from individual VARs, as detailed later, or be inferred from strike data. Moreover, it is really important to add null values in the foreland where we know structures are unrotated. Regarding the pin-element, we must select it with care and geological sense because the restoration process depends heavily on it. We should preferably select the foreland as pin-area. 3) Translation and rotation. This step is equivalent to the piecewise surface restoration method and we use the same code. Triangles are rigidly translated and rotated to minimize distances between common vertices and to fulfill the rotation constraints. An angle of free rotation is allowed, subject to the precision of data (α95), but it is better to use the rotation data as rough data. 4) Welding. After the fitting process, triangles are welded in order to have a continuous map. This step adds measurable deformation. 5) Calculus of deformation and visualization of results. We display the dilation (change in area) because it has proven to be the most representative control parameter. The anisotropy of the strain ellipse can be easily calculated from the restoration procedure. We can also display the displacement vector map, which provides relevant information as well. Moreover, we visualize the cartographic map in the restored state to facilitate the evaluation of results. This requires deforming the raster image as the restored mesh. With this purpose we have developed a C program based on openGL1, available in the companion CD with Supplementary Material. We treat the cartographic map as a texture. We triangulate this texture and then apply the proper deformation to the triangles, according to the restoration. The procedure based on the parametric approach has the same initial and final steps that the piecewise approach: 1 www.opengl.org 138 restored (β1) but not the magnitude of rotation of the second step (β2) that is proportional to the declination error. Further analysis must be done, particularly in the shear and superposed folding error, but the characterization of these errors constitutes a great help for the identification of the source of error and its later correction. 241 242 Supplementary material The work done with this thesis is an open door for further researches. Therefore, we pretend with the contents of the supplementary CD that methods and simulations can be reproduced and used as starting point for incoming studies. We detail the content of each folder. Analog models Files marzal.ts and balzes.ts are gOcad surface objects. They are the reconstruction of the analog models described in Sections 3.2.1 and 3.2.2 used to test the restoration methods. In order to restore the surface using the paleomagnetic constraint we need to select as “Default gu values” (in the Parameterizer2D window) the property with interpolated paleomagnetic vectors. Piecewise restoration The piecewise restoration code is programmed in Matlab version 6.5 and needs to be run from Matlab. We include some simulations performed in this thesis (simulation_sanMarzal.m, simulation_balzes.m and simulation_mapView.m), the fold restPmag that contains the Matlab functions and the fold ExampleData with the files of initial surfaces to restore. To run the simulations we need to include this fold in the workpath. We detail the data and functions organization and flow: Initial data: 1) Triangulated initial surface: nodes (X,Y,Z) and vertices of triangles (tri). 2) Polarity of triangles of initial surface (normal polarity: pol = 0). 3) Pin-element (starting point or line: pin = 1). 243 4) Paleomagnetism (no data value: pmag0 (dec,inc) = Inf, Inf) and its confidence angle (by default: where paleomagnetism is defined alfa95 = 0 and where is undefined alfa95 = 360). Reference unrotated data. Restoration process: 1) Optional; extrapolation of scattered paleomagnetism: [pmag_xyz,pmag_decInc]=extrapolatePmagB(tri,X,Y,Z,tri2,X2,Y2,pmag0,ref,pol); pmagExtrap=mean(pmag_decInc); decPmag=pmagExtrap(1,:,1); incPmag=pmagExtrap(1,:,2); 2) Flattening (initial paleomagnetic data is also flattened, only the resultant declination value is used in the fitting process): [tri2,X2,Y2,azimuth,dip,pmag] = flattening3(tri,X,Y,Z,pol,decPmag,incPmag); 3) Fitting (with iterating process if paleomagnetism is not hard [α95 ≠ 0] and free rotation is allowed): [tri3,X3,Y3,orderTri]= translationRotation(tri,tri2,X2,Y2,pmag,alfa95,ref,pin); n =maximum number of iterations; [X3,Y3]=iteratingB(tri,tri2,X3,Y3,alfa95,pmag,ref,n,orderTri,pin); 4) Welding: [X4,Y4] = welding(tri,tri2,X3,Y3); 5) Calculus of dilation and deformation (relation between major and minor axis and direction of major axis): d=dilation(tri2,X2,Y2,tri,X4,Y4); [rel,ang]=deformation2(tri,X4,Y4,tri3,X3,Y3); 6) Visualization of results: figure; trisurf(tri,X4,Y4,zeros(size(X4)),d); hold; axis equal; title(‘Dilation'); colormap(map1); caxis([-0.15 0.15]); view(2); triplot(tri(find(pin==1),:),X4,Y4,'r'); figure; trisurf(tri,X4,Y4,zeros(size(X4)),(rel)); hold; axis equal; title('Axis relation'); colormap(map2); caxis([1 1.3]); colorbar; 244 figure; trisurf(tri1,X4,Y4,zeros(size(X4)),(ang)); hold; axis equal; title('Maj. angle direction'); colormap(map1); caxis([0 180]); colorbar; Parametric restoration This fold contains the plugin Restoration2D for gOcad version 2009.4. This plugin is the restoration method based on the parametric approach in which we have incorporated the paleomagnetic constraint. We must set the paleomagnetic property as “Default gu value”. Map-view restoration This simple c program (restoredmap.exe) plots the cartographic map in the restored state after the palinspastic (or map-view) restoration (Section 6.1). Two files are needed to run the program: 1) a text file with the vertices of texture and position (example: vertices.txt) and a 2) bitmap file with the raster image of the cartographic map (example: carto.bmp). The program is run from the command prompt adding two arguments corresponding with these two files (example: >restoredmap vertices.txt carto.bmp). The vertices file contains the x y coordinates of texture and position. Vertices are sorted by triangles (three consecutive lines are the nodes of a triangle). Texture vertices are the points of the meshed initial map while position vertices are the points of the meshed restored map. Point range is from 0 to 1. VPD Virtual Paleomagnetic Directions program is fully described in Appendix1: Achieving the optimal information from paleomagnetic analysis. VPD.jar is the executable file whereas the source folder contains the code. We also include example files for each input format. ASN3.th and Almonacid u-channels.txt are the files used in 245 Section I.2 for the application to real data sets. We detail in the manual (VDP Manual.pdf) all routines of the program. Paleomagnetic errors These are the excel macros developed to compute the paleomagnetic errors described in Appendix 2: 1) Overlapping error Input data: P (primary vector), S (secondary), r (P/S ratio), rotation matrix (α, 0, β) Output data: Error 2) Shear error Input data: P (primary vector), bedding plane (α, β), simple shear (γ, ψ) Output data: Error 3) Superposed folding error Input data: P (primary vector), S1 (α1,β1), S2 (α2,β2) Output data: Error 246