Russian meeting in 2009

You may download some of the abstracts here or click the links below to view the others.

Interactive Visualization System: A Distributed Multiuser Web Application Using the Google Web Toolkit Framework

GPU-Based Visualization and Computational Steering of Level Sets for Detecting Subsurface Geological Features

Volumetric Visualization of Large Data Sets at LCSE in University of Minnesota

Tree-based Methods and Fast Multipole Methods on GPUs

Modeling Tsunamis with Graphics Accelerated Hardware (GPU) and Radial Basis Functions (RBF)

Some Issues on Rheological Properties of Minerals Pertinent to Global Geodynamics

Проблема количественной оценки цикла суперконтинента в фанерозое.

Mantle xenoliths are non-representative of the cratonic mantle: Reconciling thermal, seismic, and petrologic evidence

A Parallel Program for Computational Modeling of Emergence of Intensive Atmospheric Vortices

Percolation of Fe-S melts in mantle silicates: Time constraints for planetary core accretion

New crustal model for Asia and Antarctic region - first steps for building global crustal model with resolution 1×1degree

Mapping of geological substance rheidity in plan and profile

Three-Dimensional Mantle Convection Code Implemented on a Graphics-Accelerator Board (GPU ) and the high Rayleigh Number Regime.

Origin of the low rigidity if the Earth's Inner Core

Evolving rock damage near large strike-slip faults

Horizontal stress field in two-dimensional mantle convection with moving continent interacting with the mantle: a comparison of isoviscous mantle case and p,T-dependent viscosity case

Continental Tectonics, Rheological Stratification and Seismicity: Insights from geodynamic numerical modeling.

Slab Geometry and the Compatibility Condition for the Kinematics of Subduction

General parameterization and resolution of geophysical inverse problems

Lattice thermal conductivity of MgO at conditions of Earth’s interior

Wave-packets and their applications in seismology

Constraints on rheology of plastoshere from space geodetic observations of deformation following large earthquakes

Use of geoinformation models for the analysis of stress-strain state of the earth's crust of the Caspian region

DEPTH ESTIMATION FOR THE BASEMENT SURFASE USING MODIFIED SPECTRAL ANALYSIS TECHNIQUE

Crustal deformation deduced from micro-gravity and geodetic data in the High Dam area, Aswan (Egypt).

Scale aspects of heat transport in the diamond anvil cell, in spectroscopic modeling, and in Earth’s mantle

Three-dimensional numerical modeling of crustal extension

Data assimilation in problems of mantle dynamics

Constraining the effects and amount of serpentinization in the oceanic lithosphere

Influence of pressure on the deformation fabrics of olivine

Some issues on rheological properties of minerals pertinent to global geodynamics

Effective Viscosity of Suspensions of Swimming Bacteria

Rheology, deformation, shear localization and strength of the lithosphere.

Classical Density Functional Theory for Hard-Sphere Liquids

Мегациклический режим термохимической конвекции земной мантии

Internal structure of icy satellites of Jupiter and Saturn and subsurface oceans

Constraining the composition and thermal state of the cratonic mantle from inversion of seismic data

The effect of climate and damage on generating plate tectonics on a planet

Comparison between mantle seismic tomography and past plate boundaries of the western Pacific

Thermodynamics of Earth Surface Instabilities

Deformation and seismicity associated with rift zones crossing passive continental margins

ON THE POSSIBILITY OF CRATER FORMATION ASSOCIATED WITH AN ASCENDING PLUME

Generation of intermediate to deep earthquakes by self-localizing thermal runaway: insights from petrological and numerical studies

Small-scale convection in the asthenosphere: A case study of foredeep basin formation

Application of adaptive wavelets methods in computational geodynamics.

3D numerical experiments assessing the necessary conditions for a plume-fed asthenosphere

The Emergence of Low Lids in SuperSized Earths INFLUENCE OF THE PRESSURE ON THE CONVECTIVE BEHAVIOUR OF EARTHLIKE PLANETS

Novel methods in computational mineral physics: insights into minerals and compounds under pressure.

Collisional granitoids formation in rheologically layered lithosphere of overthrusted structures – numerical modelin

Modeling Transform Plate Boundary in 3D: Dead Sea Transform Fault from the Red Sea to Lebanon Mountains

Pressure variations during ultra-high pressure metasomatism

Numerical modeling of Rayleigh-Tailor instability in relation to the melting and diapirism of granite magma in the crust

StGermain software environment for HPC and software maintenance.

Magnetic filed generation in Earth’s core

Modeling of stresses caused by lithospheric density inhomogeneities: delamination of Sierra-Nevada deep batholithic root

ROLE OF SEEPAGE FORCES ON FAULTING AND EARTHQUAKE TRIGGERING

Current Mantle Energetics

Calculation of site-specific carbon-isotope fractionation in pedogenic oxide minerals

On the mantle control of core convection and the geomagnetic field

Density models and structure of the Earth`s lithosphere beneath Central and South Asia constrained with variations of free mantle surface

Three dimensional numerical modeling of lithospheric dynamics coupled with mantle convection

Studying the Earth's interior based on correlations of ambient seismic noise

Seismic tomography of Eurasia

Processes of percolation in rocks

Initiation of the modern style of subduction in the Precambrian: insights from numerical experiments

Three dimensional numerical modeling of lithospheric dynamics coupled with mantle convection

Plume formation in strongly temperature-dependent viscosity fluids and the origin of the Martian hemispheric dichotomy

Fluid flow dynamic in transient growing regime of layered visco-elastic saturated sediments.

3-D spherical models of coupled mantle thermo-chemical evolution, plate tectonics, magmatism and core evolution incorporating self-consistently calculated mineralogy

Coupling of ¨Mechanical¨ and ¨Chemical¨ Compaction Using Thermodynamic Tools

EVOLUTION OF MANTLE CONVECTION WITH SELFGENERATING LITHOSPHERIC PLATES AND MOVING CONTINENTS

Discrepancy and seismic isotropy in D''

Fine upper mantle structure from S receiver functions

Precise stationary gravimetry in studies of the Earth’ inner structure and its dynamics

Modelling of reactive fluid transport in deformable porous rocks

The Fate of the Slabs Interacting with a Viscosity Hill in Mid-Mantle.

On hypothesis of hydrogen in the Martian core

Dynamics and Implications of 3-D hydrous thermal-chemical plumes in oceanic subduction zones

Crustal and deep mantle structure of North Atlantic MOR

Gravity stresses in the Central Asia crust.

Interactive Visualization System: A Distributed Multiuser Web Application Using the Google Web Toolkit Framework
Weiss, R., McLane J., Yuen, D.A., Wang, S.M.

We have created a web-based, interactive system for collaborative, multi-user visualization of large data sets that allows users in geographically disparate locations to simultaneous and collectively visualize large data sets over the Internet. This system builds upon the large-scale interactive visualization system already in place that provides high resolution visualizations to the order of 13 million pixels by Megan Damon. By leveraging asynchronous java and XML (AJAX) web development paradigms via the Google Web Toolkit (http://code.google.com/webtoolkit/), we have been able to extend interactive visualization to hand-held devices such as the OQO, Nokia N800, netbooks, and 3G enabled cell phones with javascript capable browsers via the mobile web. In the current version of our software, we have implemented a new collaborative environment, a feature that allows multiple users to focus on and work with the same data allowing discussions to take place without regard to location. We will discuss our method of implementation, specifically why we chose the GWT framework and the benefits of a flexible modular implementation. We will also discuss the Google Application Engine (http://code.google.com/appengine/) and why we have decided against using this technology until it matures further.

Graphical description of the implementation

Application UI

Top

GPU-Based Visualization and Computational Steering of Level Sets for Detecting Subsurface Geological Features
Benjamin J. Kadlec, Dept. of Computer Sciences, University of Colorado, Boulder, Colorado, 80309, U.S.A.

A unified approach is presented for the interactive exploration of seismic datasets in real-time. This method uses dynamic implicit surfaces with a real-time visualization by computing directly on a graphics processing unit (GPU). Level set methods allow implicit handling of complex topologies deformed by operations where large changes can occur without destroying the level set representation. We solve the level set equation to steer a surface according to subsurface features imaged in the data and then immediately visualize the results in 3D.

CUDA is a general purpose parallel computing architecture that allows the NVIDIA GPU to be treated like a data parallel supercomputer in order to solve many computational problems in a fraction of the time required on a CPU. CUDA is compared to the new OpenCL™ (Open Computing Language), which is designed to run on heterogeneous computing environments but does not take advantage of low-level features in NVIDIA hardware that provide significant speedups. Therefore, our technique is implemented using CUDA and results are compared to a single CPU implementation to show the benefits of using the GPU and CUDA for seismic data interpretation.

We solve a 1024^3 problem and experience over 20x speedup using the NVIDIA GeForce 295 GTX. Implementation details for mapping the problem to the GPU architecture are described as well as discussion on porting the technique to heterogeneous devices (AMD, Intel, IBM) using OpenCL. The results present a new interactive visualization system for interpreting geologic surfaces in seismic datasets such as complex geological faults, subsurface river channels, and large ore bodies. The system allows geoscientists to observe the evolution of surfaces and steer them toward features of interest , using their domain knowledge.

Top

Volumetric Visualization of Large Data Sets at LCSE in University of Minnesota
Mike Knox, Paul Woodward, and David H. Porter Laboratory of Computational Science and Engineering (LCSE) University of Minnesota, Minneapolis, Minnesota 55455 U.S.A.

The visualization system in the Laboratory for Computational Science and Engineering (LCSE) was originally developed as a system to drive the PowerWall display in 1994. This interactive display system consists of 12 projectors which deliver an overall resolution of 15 megapixels to a three panel semi-immersive rear projection wall. Over the past four years this system has been improved and developed under NSF grants with the intended purpose of providing researchers with an on-demand resource for interactive computation and simultaneous visualization.

We are currently in the process of developing the interactive computation side of this visualization system. Over this past year we purchased six IBM Cell processor Tri-Blades, these are the same blades which compose the RoadRunner supercomputer at the Los Alamos National Lab. These blades are comprised of two dual Cell processor blades connected to a single dual dual-core AMD blade. The blades are interconnected with a DDR infiniband fabric. We have also purchased 24 dual quad-core Intel Nahelem processor nodes. These nodes will each have two Telsa GPUs and 12 of the nodes will each have 24 TB of locally attached disks, all of the nodes are connected via a QDR infiniband fabric. When in use, this system can compute a simulation using a combination of the Tesla GPUs, Intel processors, and Cell processors. While the simulation is running the output data can be analyzed and visualized using the PowerWall. The unique capabilities of this system give researchers the ability to monitor the progress of their simulation, modify it if necessary, and quickly find the answers to their questions. When computing a problem of appropriate size this system can be used interactively for simulations that only take a few hours to compute. We have demonstrated this capability at the national Supercomputing conference over the past three years.

For simulations of higher resolution we have connected the above described system to the multi-thousand core machines in the Minnesota Supercomputing Institute (MSI). While a simulation is running on one of MSI’s large machines the data can be streamed out of the compute nodes, across a 20GB/s DDR Infiniband link, to the systems and disks in the LCSE. Once the data has arrived the machines in the LCSE can analyze and visualize the data. This enables researchers to interact with a very high resolution computation. Visualizing the simulation as it is running provides researchers with the ability to ensure that the simulation is proceeding as intended, and if needed, modify the computation. This saves CPU time as well as time needed to post process and visualize data after the run has completed. MSI implemented a reservation queue which lets a researcher reserve thousands of CPUs for a specific period of time. This queue allows for planning an interactive supercomputing session at a time that is convenient for the researches and their collaborators.

This system has also proven to be valuable for remote interactive supercomputing and visualization. Over the past three years a web interface has been built on top of this system which enables remote use of the system. Using this web interface a research is able to visualize their data and monitor the progress of their simulations using any javascript enabled browser. The web interface is a full featured visualization tool exposing fine camera controls and customizable color maps. This ability has been demonstrated at the national Supercomputing conference in Austin, Texas ,as well as other international conferences, like in Elm, Switzerland.

Top

Tree-based Methods and Fast Multipole Methods on GPUs
Matthew Knepley, Dept. of Computer Science, University of Chicago, Chicago, Illinois, U.S.A.

We examine the performance of the Fast Multipole Method on heterogeneous computing devices, such as a CPU attached to an Nvidia Tesla 1060C card with close to one Tflop speed. The inherent bottleneck imposed by the tree structure is ameliorated by a refactoring of the algorithm which exposes the fine-grained dependency structure. Also, common reduction operations are refactored in order to maximize throughput. These optimizations are enabled by our concise yet powerful interface for tree-structured algorithms. Examples of performance are shown for problems arising from vortex methods for fluid mechanics, and extensions to the linear Stokes problem involving creeping flow and linear elasticity are also discussed.

Top

Modeling Tsunamis with Graphics Accelerated Hardware (GPU) and Radial Basis Functions (RBF)
Jessica Schmidt (1), Cecile Piret (2), David A. Yuen (3and4), Yingchun Liu (3), Benjamin Kadlec (5), Erik Sevre (3)
1. Dept. of Applied Mathematics, Univ. of Colorado, Boulder, CO, USA
2. National Center for Atmospheric Research, Boulder, CO, USA
3. Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, MN, USA
4. Dept. of Geology and Geophysics, University of Minnesota, Minneapolis, MN, USA
5. Dept. of Computer Science, University of Colorado, Boulder, CO, USA

Time is of extreme importance when issuing tsunami warnings. For this reason, it is important to minimize the amount of time needed to issue a tsunami warning after an earthquake strikes so that those in the affected region have the most time to react. We have utilized many tools in order to facilitate this speedup so that the tsunami may be detected early on in its propagation.

First, we used a graphics accelerator board (GPU) to test the speedup of the simulation compared to the regular CPU. By putting the linear shallow water equations on the GPU for evaluation, we obtained a speedup of about 8 times greater than the CPU alone. Furthermore, we looked at a new method to solve the wave equations, called radial basis functions (RBF). RBFs are a more novel method based on a gridless and a multi-scale, adaptive framework. Therefore, we can solve the equations using less data points, a smaller grid size, and by taking larger time steps, without compromising accuracy. Thus, this RBF method working with the GPU to solve the equations, we achieved a speedup of seven times compared to the same code running strictly on the CPU.

Finally, we have also begun to look at the seismic-tsunami-atmospheric waves. In this problem, the earthquake in the ocean causes a tsunami to form, which thereby displaces the atmosphere. By the time the gravity waves have reached the ionosphere, the proxy chemical signals will have been greatly amplified and can be detected much easier by a GPS satellite system in outer space. Hickey et al [1] have examined the initial displacement of these gravity waves in a one-dimensional case. We would like to extend their work so that we may look at the gravity waves in two-dimensions, as they propagate over time, because a slight disturbance in the ocean could be detected more easily in the upper atmosphere. Thus, this method holds much promise for tsunami detection by the larger signals induced in the chemical species.

[1] Hickey, M.P., G. Schubert, and R.L. Walterscheid, The propagation of tsunami-driven gravity waves into the thermosphere and ionosphere, J. Geophys. Res., in press, 2009

Top

Some Issues on Rheological Properties of Minerals Pertinent to Global Geodynamics
Shun-ichiro Karato Yale University, Department of Geology and Geophysics

Rheological properties of minerals control the dynamics of terrestrial planets. However, experimental studies on rheological properties are challenging not only because well-controlled deformation experiments under deep planetary conditions are difficult but also because of the large number of factors that could affect the rheological properties. In this talk, I will discuss a few topics related to the rheological properties of minerals that have a close link to the global geodynamics:

1. The strength of the lithosphere: why plate tectonics on Earth?

Previous models for lithosphere strength were based on the properties of olivine. Our recent experimental studies show that orthopyroxene has much smaller creep strength under the lithospheric conditions. I will discuss the possible implications of these observations for the operation of plate tectonics on Earth but not on Venus.

2. The energy dissipation in convecting mantle: some notes on the scaling law for thermal evolution of terrestrial planets

Thermal evolution of terrestrial planets have been studied using a “parameterized convection” model where some relationship between the average temperature and energy transport by convection are used. In a conventional model, the rate of energy transport increases with temperature, but Conrad-Hager (1999, 2000) (also Korenaga, 2003) suggested that when energy dissipation of plate bending plays an important role, then such a relation might need to be modified. I will discuss the validity of such a modification and present an alternative model for energy dissipation based on the latest knowledge of rheological properties of minerals under deep mantle conditions.

Top

Проблема количественной оценки цикла суперконтинента в фанерозое.
Апарин В.П., Петроченко С.В. ИЗК, СО РАН Иркутск, e-mail palmag@crust.irk.ru

В последние двадцатилетия в геотектонике снова оживились споры о роли цикла суперконтинента (CS) в тектонике плит [1-4]. Объекты обсуждения затрагивают вопросы формирования CS, наличие возрастных и палеогеографических рубежей, блоков-участников цикла суперконтинента (CS), ранговую классификацию циклов [2,6].
При этом элементы критики моделей ЦС связаны с тем, что на финальном этапе полной сборки суперконтинента (Пангея) обнаруживаются дефекты: траектории палеомагнитных полюсов (APWp) и края материков, участвующих в аккреции не дают удовлетворительного совмещения [1-2]
В частности, задача данной работы связана с изучением скрытой периодичности во временных рядах CS в фанерозое, то есть в истории Пангеи. Возросли требования к унификации подходов и методов в работе с временными рядами (всегда использовать единую возрастную шкалу, единую систему расчетов и т.п.). В данной работе приведены некоторые результаты исследований процессов, полученные методами спектрального анализа временных рядов из материалов различных дисциплин (палеонтология, палеомагнетизм, геохимия осадочных пород, включая изотопные возрастные ряды датированных событий.) Для исследования рядов использовались два метода спектрального анализа: метод максимальной энтропии и вейвлет-анализ [ 5 ] .
Почти в каждом объекте обнаружена сложная хронологическая структура процессов, которые фиксировали суперпозицию стационарных компонент и детерминированного тренда [5]. Чтобы детально изучить скрытую периодичность, извлекались тренды. В результате чего удалось установить скрытую периодичность палеонтологической летописи, которая фиксирует периодичность развития органического мира океана. В итоге в спектрах макроэволюции биоты установлено: а) доминирование стационарных компонент циклов с длительностью периодов 364 и 400 млн. лет. б) распределение экстремумов на возрастных зависимостях доминирующих компонент образует рисунок, устойчивый на протяжении фанерозоя в) обнаружен также незначительный вклад высокочастотных компонент в энергетику спектров временных рядов морской фауны [4,7].
Количественные оценки истории чехла Восточно-Европейской платформы (ВЕП) представляет эталон изученности временных рядов осадочного чехла всех древних платформ планеты. Погрешность оценок не выходит за рамки 10%. [2].. Процесс образования осадочного чехла ВЕП в фанерозое (спектральный анализ) подчиняется периодичности в области низких частот, где доминируют циклы с периодами 200 и 220млн. лет для вариации площади морей ВЕП. Для других платформ эти цифры близки к 300 млн. лет. Геохимические показатели чехла ВЕП для глин и алевро-песчаников подчиняются временным вариациям тем же, что общие палеогеографические параметры чехла: так процентное содержание С орг., СО2, СаО, во временных рядах осадков чехла ВЕП также подчиняются циклу 200 млн. лет. Процентное содержание Аl2О3 и К2О в алевролитах и песчаниках изменяется с периодичностью 400 и 364 млн. лет. Совмещение полученных гармоник с исходными временными рядами показывает сопряженность главных экстремумов на возрастной шкале для многих изученных временных зависимостей осадочного чехла ВЕП. В частности вариация и форма движения на цилиндрической проекции показывает связь движения плит с накоплением осадочного материала, а петли APWp характеризуют [4] амплитуду тектонических деформаций.
Сходство динамики фанерозойского осадконакопления на платформах ВЕП и САП отмечено в прошлом веке [2]. Главный фактор этого сходства образование Еврамерики в финале каледонского и герцинского тектонических циклов и закрытие Уральского палеоокеана [6,7] Наложение временных рядов вариаций площади моря древних платформ [2,6] показывает сопряженность позиций экстремумов ВЕП и САП, не смотря на различие в длительности циклов. Стратиграфические границы блоков фиксируются от перерыва до перерыва в геологическом разрезе. Динамика возрастных вариаций числа блоков на геохронологической шкале отражает роль тектогенеза в истории чехла САП. По данным спектрального анализа временных рядов распределения числа блоков в фанерозое доминирует цикл с продолжительностью 400 млн. лет., экстремумы которого на возрастной шкале близки циклу морского осадконакопления и рубежам APWp Северной Америки в фанерозое.
В итоге, даже простое перечисление истории природных объектов в составе CS (литосферных блоков) показывает, что преимущество имеют горизонтальные перемещения блоков за счет конвекции в мантии.

Литература
1. Condie K.C. The supercontinental cycle: are two pattens of cyclicity // Jorn. of African Earth Sci.-2002.- v, 53.-p179-183.
2. Добрецов Н.Л., Кирдяшкин А. Н., Кирдяшкин А.А. Глубинная геодинамика // Новосибирск,- изд-во СО РАН.- филиал «ГЕО», 2001,- с.408.
3. Апарин В.П., Альдуганов Д.Н. Периодичность колебаний уровня океана и распределения запасов нефти в фанерозое. //ДАН РАН, 2001, т. 376, №5, с.684-688.
4. Апарин В.П., Петроченко С.В., Кузьмин А.В., Кашкин В.Б. Низкочастотные компоненты в вариациях разнообразия морской биоты и уровня мирового океана.// ДАН РАН, 2004.т.397, с. 108-112.
5. Мала С. Вейвлеты в обработке сигналов. //М., «Мир»,2005, с.452
6. Трубинин В.П., Рыков В.В. Проблемы глобальной геодинамики. М. «ГЕОС»,2000,-с.7-28.
7. Лисицин А.П. Лавинная седиментация и перерывы осадконакопления в морях и океанах //М. Наука, 1988, - с.309.
8. Апарин В.П.,Петроченко С.В. Длиннопериодные колебания уровня океана и цикл суперконтинента в фанерозое. //ДАН РАН, 2005, т.405 №4 с.539-543
9. Апарин В.П., Сарычев В.Т. Связь длиннопериодных вариаций частоты геомагнитных инверсий с осадочными циклами в фанерозое. // Геомагнетизм и аэрономия 2000, т.40, с.141-145
10. Апарин В.П., Кузьмин А.В. Оценка факторов периодичности карбонатонакопления в фанерозое. // ДАН РАН 2006, т. 408, №5

Top

Mantle xenoliths are non-representative of the cratonic mantle: Reconciling thermal, seismic, and petrologic evidence
Irina M. Artemieva University of Copenhagen, Denmark, Irina@geo.ku.dk; www.lithosphere.info

I use three independent approaches to argue that xenolith data on densities and seismic velocities in the Archean-Proterozoic lithospheric mantle maybe non-representative of the "intact" (=unsampled by xenoliths) cratonic mantle.

The first example (Artemieva, 2007) is based on buoyancy modeling for the East European Craton (EEC), that lacks surface relief but has huge amplitudes of topography at the top of the basement (20+ km), at the crustal base (ca. 30 km), and at the lithosphere-asthenosphere boundary (200+ km). I examine the relative contributions of the crust, the subcrustal lithosphere, and the dynamic support of the sublithospheric mantle (SLM) to maintain surface topography, using regional seismic data on the crustal structure and thermal data. The results of buoyancy modeling indicate either a smaller density deficit (ca. 0.9 per cent) in the SLM of the Archean- Paleoproterozoic parts of the EEC than predicted by global data on mantle xenoliths (ca. 1.5 per cent) or the presence of a strong convective downwelling in the mantle beneath the craton.

The second example (Artemieva, 2009) is based on geophysical constraints on global-scale compositional variations in the SLM. I calculate seismic velocity variations of a non-thermal origin in the SLM of the continents (using global Vs surface-wave and body-wave seismic tomography data and lithospheric temperatures). In agreement with xenolith data, the results show strong positive velocity anomalies of non-thermal origin (attributed to mantle depletion) for all of the cratons; their amplitude varies laterally and decreases with depth. SLM of the cratonic regions where kimberlite magmas erupted show only weakly positive compositional velocity anomalies as compared to the adjacent "intact" cratonic mantle, implying that xenolith data maybe non-representative of the cratonic mantle.

The third example (Kaban et al., 2003) is based on the results of global gravity modeling in which the effect of spatially differential thermal expansion has been eliminated from mantle residual gravity anomalies. These results indicate a large scatter of density deficit in the cratonic lithosphere that is not correlated with crustal differentiation ages, in contrast to global xenolith data.

I conclude that variations in mantle densities and velocities constrained by xenolith data and parameterized solely in terms of iron-depletion are non-representative of "intact" cratonic mantle. The discrepancy between major-element petrologic models and geophysical observations stems from at least two sources: (a) metasomatic modification of the cratonic lithosphere prior/during kimberlite magmatism; (b) underestimate of the effect of orthopyroxene on densities and velocities.

References:
Artemieva I.M., 2009. The continental lithosphere: Reconciling thermal, seismic, and petrologic data. Lithos, 109, 23-46.
Artemieva I.M., 2007. Dynamic topography of the East European Craton: Shedding light upon the lithospheric structure, composition and mantle dynamics. Global Planet. Change, v. 58, 411-434.
Kaban M.K., Schwintzer P., Artemieva I.M., and Mooney W.D., 2003. Density of continental roots: compositional and thermal effects. Earth Planet. Sci. Lett., 209, 53-69.

Top

A Parallel Program for Computational Modeling of Emergence of Intensive Atmospheric Vortices
Babkova B.V. Institute for System Programming of RAS,Moscow, barbara@ispras.ru,
Gubar A.U. Schmidt Institute for Physics of the Earth of RAS, Moscow, parkAG@yandex.ru

A scalable parallel program for solution of nonlinear system of equations, modeling processes and start conditions of Intensive Atmospheric Vortexes (IAV) in three-dimensional compressible atmosphere is discussed. The system of equations vas obtained in [1,2] and is a strongly nonlinear system of the mixed type. The program was developed in the Institute for System Programming of RAS in collaboration with the Institute for Physics of the Earth of RAS using ParJava[3] Environment and JavaMPI library. An explicit conditionally stable scheme of the second order of time and space accuracy was used, its stability criterion was close to that of the explicit MacCormack scheme [4].

When a parallel program is being developed or modified, it is necessary to make sure not only its correctness but also its scalability. However, the analysis of dynamic properties of a program (profiles, traces, slices, etc.) allowing to state its scalability, usually is coupled with necessity of numerous running of the program, which is not completely debugged yet, on target computing system (high-performance cluster). The tools of ParJava Environment allowed to perform the most part of the analysis of a parallel program using an instrumental computer. This reduced the duration and computational burden needed for development and debugging of the code.

The program can be divided in two parts: loading/initialization of data and the main loop. Saving of the results occurs during the main loop execution. The program’s input data (the physical parameters of the model and housekeeping data) are stored in the separate file.

After initialization the program saves zero data layer and starts the main loop. Layer() function, computing the values of the main arrays, is called four times during each iteration

All loops were explored using distance test to reveal the possibility of their parallelizing. The absence of dependencies provided to divide arrays into blocks and distribute the blocks among processors of a cluster. Since a difference scheme is used, the blocks must overlap. The regions of overlapping (shadow faces) during the computation must be passed from the processor computing them to the processors using them. For the three-point scheme the width of a shadow face equals to one space layer. During the current iteration only data computed during previous iterations are used. It allows to refresh shadow faces only once for each level cutting overheads. As was shown by analysis two-dimensional partition of the array is more efficient than one-dimensional one, so two-dimensional partition was used.

Significant results are stored by each process on a local hard disk in binary arrays. Merging of the data stored during program execution is performed after the end of computation. Passing of such volumes of data during execution of the program would cause strong mistiming on each step and would considerably augment the total execution time.

Since the program of tornado simulation is a large-scale program both in running time and in data, the modified checkpoint subsystem was used/ It allowed to decrease the volume of the task context stored by 25%. For reason of saving of disk space we decided to resample down the arrays with visualization data. Being unessential to the quality of visualization it allowed to decrease in 8 times the space needed for each set of visualization data.

Since the program of tornado simulation is a large-scale program both in running time and in data, the modified checkpoint subsystem was used/ It allowed to decrease the volume of the task context stored by 25%. For reason of saving of disk space we decided to resample down the arrays with visualization data. Being unessential to the quality of visualization it allowed to decrease in 8 times the space needed for each set of visualization data.

Fig. 1 Speed-up graphics.

Presented results of computations were used in research of the process of tornado origination and demonstrated adequacy of the model being used and workability of ParJava environment for development of such applications as well.

LITERATURE: 1. Arsenyev S.A., Gubar A.Yu., Nikolaevskiy V.N. Self-Organization of Tornado and Hurricanes in Atmospheric Currents with Meso-Scale Eddies. // Doclady Earh Science., 2004, т.396, №4, с.541-546
2. Nikolaevskiy V.N. Vortexes and waves. M: Mir, pp.266 – 335, 1984..
3. Victor Ivannikov, Serguei Gaissaryan, Arutyun Avetisyan, Vartan Padaryan. Improving properties of a parallel program in ParJava Environment // The 10th EuroPVM/MPI conference. LNCS 2840. Sept. 2003, Venice. pp. 491-494.
4. Fletcher C.A.J. Computational Techniques for Fluid Dynamics. Springer-Verlag Berlin Heildeberg. 1991.

Top

Percolation of Fe-S melts in mantle silicates: Time constraints for planetary core accretion
N. Bagdassarov Institute of Geosciences, University Frankfurt, Germany; nickbagd@geophysik.uni-frankfurt.de

The mechanism which segregates molten Fe-S into metallic cores of planetary bodies is still not fully understood. Due to the high interfacial energy and wetting angle between Fe-S melts and silicate mantle minerals, continuous percolative flow of such melts cannot be efficient for core segregation in planetary bodies. A series of percolation experiments has thus been realised on partially molten fertile garnet peridotite, employing a centrifuging piston-cylinder. The high temperature garnet-peridotite with Mg# ~ 0.90 is composed of 60 vol% olivine, 15 vol% orthopyroxene, 6 vol% clinopyroxene and 19 vol% garnet has been used as a silicate matrix. Peridotite powders with 100-200 or 20-30 µm grain size were mixed with 5-30 vol% Fe-S of eutectic composition. The aggregates were centrifuged at 500-700 g at temperatures below and above the melting point of peridotite. The centrifuge experiments revealed a negligible percolation of Fe70S30 melts through the unmolten peridotite matrix. Only at T>1260° C, i.e. above the solidus of the peridotite, and starting with 5 vol% of Fe70S30 the vertical melt gradient achieved 1-2 vol%/mm. In samples with 15 vol% Fe70S30 the vertical separation achieved 2-2.5 vol%/mm after 10h of centrifuging at 500 g. An increase in the degree of partial silicate melting in the peridotite leads to an increase of the Fe70S30 separation rate from the peridotite matrix. Fe-S contents > 10 vol% cause an increase of the Fe70S30 melt droplet size and of the effective percolation velocity of Fe70S30 melt, in agreement with the results of Yoshino and Watson (2005). The percolation threshold dividing the fast (>10 cm per year) and slow percolations (<1 mm per year) of Fe70S30 melt is around 14-15 vol% of Fe70S30. The experimentally determined permeabilities of Fe70S30 melt in unmolten peridotite with 7-10 vol% of Fe70S30 melt are 10-17-10-18 m2, which is 1-2 orders of magnitude lower than the values calculated by Roberts et al. (2007) from static experiments. The presence of silicate melt increases the segregation velocity of Fe70S30 in partially molten peridotite by more one order of magnitude with respect to unmolten peridotite matrix. This could provide an effective segregation of Fe70S30 melt in a planetary mantle down to 2.5-3.0 vol% of residual Fe-S melt. The extremely slow percolation of Fe70S30 melt in the absence of partial silicate melting precludes a scenario of metallic core formation via percolation before temperatures allow for partial melting of mantle silicates in planetary bodies.

The connectivity of FeS-melts in olivine and in a fertile peridotite matrix has been also addressed through in-situ electric impedance spectroscopy (IS) measurements at 1 GPa. A first series of experiments used sintered powder samples of a fertile peridotite xenolith mixed with 5-15 vol% Fe70S30 of eutectic composition. The sheared high-T garnet peridotite with Mg# ~ 0.90 is composed of 60 vol% olivine, 15% orthopyroxene, 5.3 % clinopyroxene and 19% garnet, the powder grain size was 20-30 µm, similar to the one employed by Yoshino et al. (2004). For a second series, San Carlos olivine aggregates were used as solid matrix and 10-20 vol.% of eutectic Fe70S30 were added. For these, the average grain size was 3 µm, much smaller than in the experiments by Yoshino et al. (2003). The powder mixtures of peridotite + Fe70S30 and olivine aggregates + Fe70S30 were first annealed for 2-5 days in a conventional piston-cylinder at 1 GPa, 950-970°C. The electrical conductivity of samples was then measured using the impedance spectroscopy method in BN-graphite-CaF2 pressure cell with concentric cylindrical electrodes made from Mo- and Re-foil (the assembly corresponds to an oxygen fugacity close to the IW buffer). The results indicate that up to 15 vol% of Fe70S30 the melt phase does not built a stable interconnected network in a peridotite matrix, as was recently indicated by Walte et al. (2007). The percolation threshold for a stable FeS network in olivine matrix lies at 17.5 vol%, much higher the 6 vol% found by Yoshino et al. (2003). Our result is in line with the high dihedral angles of typically 90-100o for Fe-S melts in mantle materials. The higher interconnectivity threshold of this study, as compared to previous studies (Yoshino et al., 2003, 2004; Roberts et al., 2007) is a result of our smaller starting grain sizes (for olivine) in combination with much longer run durations. Both these experimental conditions result in enhanced grain growth and thus to a higher degree of textural equilibration, leading to the occurrence of the time depending pinging off of Fe-S melt films in our texturally more mature experiments.

Fig. 1. Measured segregated velocities of silicate melt and Fe-S melt in fertile peridotite. Upper panel: The scarce data of the segregation velocity as a function of silicate melt content are attributed to the differing content of Fe-S melt in samples. Dashed line stands for the correlation of silicate melt segregation velocity vs silicate melt fraction : Vsilicate melt(mm/h)=2∙10-6∙e0.33∙ Solid line represents the correlation between segregation velocity of Fe-S and silicate melt fraction: VFe-S(mm/h)=10-5∙e0.26∙. Lower panel: The correlation between segregation velocity and Fe-S content is more obvious. Dashed line stands for the correlation of silicate melt segregation velocity vs Fe-S melt fraction : Vsilicate melt(mm/h)=2∙10-5∙e0.27∙. . Solid line represents the correlation between segregation velocity of Fe-S and silicate melt fraction:VFe-S(mm/h)=5∙10-6∙e0.34∙ Horizontal arrows indicate a typical migration velocity of Fe-S in order to form metallic core by percolation process. Segregation velocities for Vesta and Mars are reduced to the Earth’ conditions.

Top

New crustal model for Asia and Antarctic region - first steps for building global crustal model with resolution 1×1degree
A.A. Baranov Schmidt Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russian Federation (baranov@ifz.ru)

Many people in different parts of earth sciences use now global CRUST2.0 model with resolution 2x2 degree. But for present days this model stays rather poor and has not enough resolution for accurate calculations. When we understand this we decide to improve this model using last years seismic data which stay available for different continents and try to build new global crustal model with resolution 1x1 degree.

We decide to begin from Asia as a big and complex continent in tectonics. The Southern and Central Asia is tectonically complex region with great collision between Asian and Indian plates and its evolution is strongly related to the active subduction along the Pacific border. Model AsCRUST-08 (Baranov,2008) of Central and Southern Asia with resolution 1x1 degree sufficiently improved CRUST2.0 in several regions and we built integrated model of the crust for Asia region. Also we add several regions in North Eurasia as Mongolia, Kazahstan and others. For such regions as Red and Dead sea, Northern China, Southern India we built regional maps with more detailed resolution. It was used data of deep seismic reflection, refraction and receiver functions studies from published papers. The existing data were verified and crosschecked. As the first result, we demonstrate a new Moho map for the region. (fig.1) The complex crustal model consists of three layers: upper, middle and lower crust. Besides depth to the boundaries, we provide average P-wave velocities and densities in the upper, middle and lower parts of the crystalline crust. Limits for Vp velocities are: for upper crust 5.5-6.2 km/s, for middle 6.0-6.6km/s, for lower crust 6.6-7.5km/s. We recalculated seismic P velocity data to density in crustal layers using rheology properties and geology data.

Fig 1. Depth to Moho for Central and Southern Asia.

As the second we selected Antactica as a child of Gondwanaland to compare crustal models for parts of Laurasia (Asia region ) and Antarctica Different tectonic units cover the Antarctic territory: platform, orogen and depression structures. This structural variability is reflected both in thickness and physical properties of the crust. Previous crustal model (CRUST 2.0.Bassin et al. 2000) for Antarctica don’t meet present-day requirements. A lot of new seismic data and regional compilations became available during last several years. We used data of deep seismic reflection, refraction and receiver functions studies as well as existing regional models (e.g. for Maud Land region, Hoffmann et al., 2003) from published papers and integrate them in a new model at a uniform grid with resolution of 1x1 degree. A new digital 3D model for the crust of Western and Eastern Antarctica and surroundings have been built. The existing data were verified and crosschecked. We present a suite of crustal models within the main tectonic units:West Antarctica rift system (WARS), the Transantarctic Mountains (TAMs), and East Antarctica (EA). As for Asia we demonstrate a new Moho map for the region. (fig.2)

Fig 2. Depth to Moho for Antarctica.

The new map demonstrates the large differences with previous models. It turns out that many regions are more heterogeneous than it was demonstrated by the previous compilations. The crustal model comprises 3 layers of crystalline crust. For each of the three basic layers the thickness and the P-wave seismic velocity (Vp) are displayed. The West Antarctic rift system is one of the largest zones of continental extension on the Earth. The seismic data show a thin extended continental crust. Crustal thickness of WARS is variable from 21 km in the Bentley subglacial trench, to 32 km in the southern flank of the Marie Byrd Land. Transantarctic Mountains: 4000 km long, peaks 4 km above Sea Level, 200-300 km wide. TAMs are characterized by the rather strong variations of the Moho depth (28-40 km). Further inland, beneath the TAM, the estimated Moho depths range from 30–33 km (30 km from the coast) to 36–40 km (85 km from the coast), deepening away from the coast beneath the TAM. Crustal thickness of EA is variable from 31 km in the Wilkes basin, to 50 km in the Western Maud Land Region.

Conclusions:
Moho map and the velocity structure of the crust are much more heterogeneous than in previous maps CRUST 2.0 (Bassin et al., 2000), and CRUST 5.1 (Mooney et al., 1998). Our model offers a starting point for numerical modeling of deep structures by allowing correction for crustal effects beforehand and to resolve trade-off with mantle heterogeneities. So we do first steps for building global crustal model with resolution 1 x 1 degree. This model will be used as a starting point in the gravity modeling of the lithosphere and mantle structure. In future we need to build such regional models for South America, Africa, Australia and North Asia.

[1] A. Baranov, Complex 3D crustal model of Asia region, Geophysical Research Abstracts, Vol. 11, EGU2009-984-5, 2009, http://meetings.copernicus.org/egu2009/
[2] A. Baranov, AntarcticCRUST-08: New crustal model of Antarctica region based on seismic data - next step for building global crustal model with resolution of 1 x 1 degree, Geophysical Research Abstracts, Vol. 11, EGU2009-1792-2, 2009, http://meetings.copernicus.org/egu2009/
[3] Bassin et al., The Current Limits of Resolution for Surface Wave Tomography in North America, // EOS Trans AGU. 2000. 81(48), Fall Meet. Suppl., Abstract F897, http://mahi.ucsd.edu/Gabi/rem.html
[4] M. Hoffmann, A. Eckstaller, W. Jokat, H. Miller, 2003. Development of a 3-D Crustal Model in the Western Dronning Maud Land Region, Antarctica, from the Interpretation of different geophysical data sets.

Top

Mapping of geological substance rheidity in plan and profile
Baranov I.P. Institute of Biological Instrument-making, RAS, Pushchino, CEK-MO@rambler.ru

Many substances are known to be able to display flow in a solid state. Under strains for several days ice begins to deport itself likes liquid. As a result of deformation forces for more than ten years salt becomes similar to flow and under action of these forces during more than 100 thousand years granite gets fluidity as a feature. But the best examples of rheid flow are crystalline shales and gneisses which have been originally deformed in the inner parts of the crust. Within a sufficient period of time all substances demonstrate flow and undergo plastic deformation in a solid state. This type of deformation is called rheidity.

Carey (by King [1]) has found that ancient strata can be contorted and extruded and resemble flow identical in every respect to flow of viscous fluid. King [1] thinks that “time is a decisive factor, and since geological time is a period of high order, all the Earth on the whole may be considered as ultimately demonstrating properties of fluid body. Under a short time strain as seismic impact and breaks, its big tectonic blocks may display itself as solid bodies”. At the same time a group of scientists headed by Prof. I.N.Stepanov in Pushchino Research Center succeeded in mapping most exactly and credibly geological substance flow in plan. They were the first who suggested a model of formalized, mathematical mapping of surface relief of geological horizons as flow-oriented lithodynamic structures. The basis of this model are isolines of zero curvature which reflects geological space curvature in earth gravity field. As a result of transformation of isolines and isohypses of topographic and structural maps, respectively, relief appears to be two systems of treelike structures – positive and negative, which differ from each other in the rate of flow of lithologic masses from which these structures are composed of. Unlike negative structures positive flow-oriented ones have greater masses and are exposed to considerable strains in gravitation field and consequently their rate is faster. Flows like these are formed as a result of mass distribution in the Earth gravity field both on a surface and in all geological horizons of the Earth’s crust. Sediments, lithodynamic flows which lie hidden at a depth of tens, hundreds and thousands meters, continue their subsurface flow extended in time over thousand and million years depending on the lithologic composition. Natural scales provide an idea of the age of one or another flow structures. The more ancient and deep structures are best depicted in maps of small scale where negligible unevenness of the Earth surface are not seen. The bigger is natural scale, the more detailed the younger rheid structures, in the formation of which deeply involves surface water, are mapped. These cartographic models are named by Academician V.P.Volobuev as “relief’s plasticity maps”.

In parallel with this model the author has used geometric analysis of curvature of geological space by seismic profile. Due to this approach it became possible to detect seismodynamic flows by profile curvature. These seismodynamic flows are formed in the author’s opinion as a result of differentiation of geological substance in profile. It implies that every profile is a fixed picture of flows of geological substances of different density towards a surface of the Earth at particular instants of time.

Methods alluded may be used in many branches of geology including oil geology. For instance, inhomogeneity of geological space in plan at a depth of more than 3 km and also in profile is shown in Figure using an oil field in Krasnodar Territory as an example. Two efficient wells shown as black spots in plan are located within a large positive lithodynamic structure being not of necessity in the center of hill crown. Most likely within planned lithodynamic structures as a result of high surface strain favorable conditions are generated for keeping basic collectors. At the same time by vertical seismodynamic flows as by channels from the region of generation fluids of hydrocarbons reach pool rock keeping in oil traps depicted in plan. These vertical flows have to continue their flow. As shown in the right side of the Figure of seismological model, flows which keep direction but have discontinuity in their flow are no longer perspective.

Figure. A complex model of geometric analysis of geological space in plan and profile using a oil field in Krasnodar Territory as an example

Despite the fact that most scholars hold the traditional notion that the Earth is static and moved only by blocks like slices of puff pastry, the above methods made it possible to confirm rheid properties of the strata of the Earth’s crust and represent them as discrete, lithodynamic flows, flow structures and flow systems. Therefore, all land surface and Earth crust are composed of treelike flow structures.

Reference
1. King L. «Morphology of the Earth» (Study and synthesis of data of the relief of the Earth). «Progress», Moscow, 1967

Top

Three-Dimensional Mantle Convection Code Implemented on a Graphics-Accelerator Board (GPU ) and the high Rayleigh Number Regime.
Gregory A. Barnett, Jr., Dept of Applied Mathematics, University of Colorado, Boulder, Colorado, U.S.A.
Grady B. Wright, Dept. of Mathematics, Boise State University, Boise, Idaho, U.S.A.
David A. Yuen, Dept. of Geology and Geophysics and Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota, U.S.A.

The last decade has seen the strong influence exerted by the gaming industry on high-performance scientific computing since the landmark year of 2003 when GPU’s speed of floating point operations surpassed that of CPU. From that time on, the amazing progress of GPU has accelerated greatly and this progress is very astounding because of the development of faster components with many more flow processing units with larger memories, now up to the 4 Gbyte range in total . There is at least a factor of 10 between the raw working speeds of GPU and CPU because there are many more processing elements ( cores ), up to 1000 cores in the case of the Tesla, in a GPU . Now with the Nvidia Tesla box it is feasible to have a potential computing power of around four Teraflops in your office environment for around 5000 dollars .

Up to now, there have already been many codes for hyperbolic problems in astrophysics, molecular dynamics, quantum chemistry, seismology, tsunami waves ported over to GPU systems. There still remain the challenges for non-linear parabolic and non-linear elliptic systems, because of the bandwidth limitations associated with the sparse matrix system of the elliptic operator, which is endemic for both CPU and GPU hardware. In this talk we will focus our attention on mantle convection, where there are both elliptic component in the momentum equation and parabolic component in the temperature equation. The mantle convection equations we have taken come from the paper by Tine Larsen et al. ( 1997), where a constant viscosity has been assumed and we have used two scalar potential functions, omega and phi to solve the momentum equations by casting them into two coupled 3-D Poisson equations in which the temperature perturbation is the forcing term of the omega potential. The poloidal velocity is derived by taking the double curl operator of the potential phi, which is obtained by solving a Poisson equation with the omega potential as the heterogeneous forcing term on the right-hand side.

All spatial derivatives are approximated by fourthorder compact finite-difference operators, with coefficients derived using Taylor series expansions. For example, to approximate the first derivative with respect to z at each node in an mby-n-by-p 3-D array f, one solves the following system of equations for f'.

h*[1/4*f'(i,j,k-1) + f'(i,j,k) + 1/4*f'(i,j,k+1)] = -3/4*f(i,j,k-1) + 3/4*f(i,j,k+1), for i from 1 to m, j from 1 to n, k from 2 to p-1,

h*[-6*f'(i,j,1) -18*f'(i,j,2)] = 17*f(i,j,1) - 6*f(i,j,2) - 9*f(i,j,3) + f(i,j,4), h*[18*f'(i,j,p-1) + 6*f'(i,j,p)] = f(i,j,p-3) - 9*f(i,j,p-2) - 9*f(i,j,p-1) + 17*f(i,j,p), for i from 1 to m, j from 1 to n.

In all cases the associated linear system is tridiagonal, and is solved with the vectorized Thomas Algorithm, requiring O(m*n*p) operations.

A Fourier Method is used to solve Poisson's equation (the momentum equation) by first looking at the linear system which results from naturally re-ordering the array into one long column-vector. More specifically, since the boundary conditions in the z direction are Dirichlet in nature, a fast sine transform is applied to the input array in the z direction to decouple the (m*n*p)-by-(m*n*p) linear system into p (m*n)-by-(m*n) linear systems. Likewise, because the remaining boundary conditions are Neumann in form because of the reflective nature of both the horizontal heat flux and the shear stress components, a fast cosine transform is applied to the columns which further breaks each (m*n)-by(m*n) system into n m-by-m systems, each of which is tri-diagonal. The boundary conditions are also fourth-order correct in all directions.

The energy equation is advanced through time by a third order explicit O.D.E. Runge-Kutta method. That is, after solving the momentum equations and approximating the spatial derivatives, the energy equations reduces to the following ordinary differential equation by the method of lines, where T is a column vector of the order of 10**6 to 10**8 long ,representing the discretized temperature field at each of the 3-D grid point

dT/dt = A(T) - B(v,T),

where A and B are functions which approximate the Laplacian of T, and v dot grad(T), respectively. Here v is the velocity vector given by the velocity potential phi in Larsen et al. (1997).We show in the figure below the flow chart of the numerical scheme where the time-stepping is illustrated.

In the connection with GPU, we have availed ourselves of the Cuda FFT ( Fast Fourier Transform) routines on the the Nvidia GPU, which facilitates the rapid solution the matrix equations associated with the Poisson equations on the GPU. We have implemented this 4th order correct algorithm on a Nvidia GTX-295 GPU by using Matlab routines, Accelereyes package for conversion to C language and finally a CUDA compiler by Nvidia. We found a speed-up factor of about 10. The computational time for a grid-size of 2 million grid points /time step is around 0.1 sec for one GPU. We will report the first GPU derived mantle convection results for 3-D solutions with Rayleigh number up to 100 million . An aspect-ratio of 2x2x1 has been employed throughout with unity being the depth.

Larsen, T.B., Yuen, D.A., Moser, J. and B. Fornberg, A high-order finite-difference method applied to large Rayleigh number mantle convection, Geophysical and Astrophysical Fluid Dynamics, 84, 57-77, 1997.

Top

Origin of the low rigidity if the Earth's Inner Core
Anatoly B. Belonoshko Applied Materials Physics, Dept of Materials Sci & Engn, The Royal Inst of Technology, 100 44 Stockholm, Sweden e-mail: anatoly@fysik.uu.se

Earth's solid-iron inner core has a low rigidity that manifests itself in the anomalously low velocities of shear waves as compared to shear wave velocities measured in iron alloys. Normally, when estimating the elastic properties of the core, one calculates an average over different orientation of the ideal iron/iron alloy crystals. This approach does not take into account the defects that are likely to be abundant at high temperatures relevant for the inner core conditions. By using molecular dynamics simulations, I show that, if defects are considered, the calculated shear modulus and shear wave velocity decrease dramatically as compared to those estimates obtained from the average single-crystal values. To demonstrate that, I first grow a crystal of iron where the equilibrium number of defects is present. This crystal consists of 16 millions atoms simulated with embedded atom model. Then this crystal is deformed and shear modulus is computed as the elastic response to deformation. The response is in agreement with the low rigidity of the inner core. Thus, the low shear wave velocity the inner core explained: it is due to the defects of the ideal lattice.

Top

Evolving rock damage near large strike-slip faults
Yehuda Ben-Zion Department of Earth Sciences, University of Southern California, Los Angeles, CA, 90089-0740, USA

Earthquake ruptures and seismic radiation generate near-fault stresses larger than the elastic limit. This produces increasing crack density (rock damage), leading to reduction of seismic velocities and increasing attenuation. During the interseismic periods, regaining of cohesion produces rock healing. Since the damage generation is resisted by normal stress, and healing is enhanced by it, the rock damage is significant in the shallow crust. Analysis of seismic data recorded near the North Anatolian fault shows clear regional reduction of seismic velocities in the top 100-500 meters during the time of the 1999 Mw7.1 Duzce earthquake (Peng and Ben-Zion, 2006). Spectral ratios of seismograms recorded at fault-zone and off-fault sites show changes of spectral curves right after the Duzce earthquake consistent with a reduction of the S-wave velocity at the fault-zone site of 20-50 per cent (Wu et al., 2009). The co-seismic changes are followed by a logarithmic recovery that is very pronounced in the first day and continues with appreciable amplitude in the subsequent 3 months. The observations are compatible with results of laboratory experiments with granular material and rocks (e.g., Johnson and Jia, 2005; Dieterich and Kilgore, 1996), and they can be simulated with a nonlinear continuum damage rheology (Lyakhovsky et al., 2009). The damage generation and healing over many earthquake cycles produces a flower-type fault zone structure, with significant shallow damage that decreases in amplitude and width with depth. The pronounced shallow damage has a hierarchical structure, consisting of ~1-3 km broad zone with enhanced seismic scattering and fault-parallel anisotropy (e.g., Peng and Ben-Zion, 2004, 2005), and more intense ~100-200 m wide damage zone that can trap seismic waves (e.g., Li et al., 1994; Ben-Zion et al., 2003) and may be manifested near the surface by a belt of pulverized rocks (e.g., Dor et al., 2006).

References
Ben-Zion, Y., Z. Peng, D. Okaya, L. Seeber, J. G. Armbruster, N. Ozer, A. J. Michael, S. Baris and M. Aktar, A shallow fault zone structure illuminated by trapped waves in the Karadere-Duzce branch of the North Anatolian Fault, western Turkey, Geophys. J. Int., 152, 699-717, 2003.
Dieterich, J.H. & Kilgore, B.D., 1996. Imaging surface contacts: power law contact distributions and contact stresses in quarts, calcite, glass, and acrylic plastic, Tectonophysics, 256, 219-239.
Dor O., Y. Ben-Zion, T. K. Rockwell and J. Brune, Pulverized Rocks in the Mojave section of the San Andreas Fault Zone, Earth Planet. Sci. Lett., 245, 642-654, doi:10.1016/j.epsl.2006.03.034, 2006.
Johnson, P.A. & Jia, X., 2005. Nonlinear dynamics, granular media and dynamic earthquake triggering, Nature, 473, 871–874.

Top

Horizontal stress field in two-dimensional mantle convection with moving continent interacting with the mantle: a comparison of isoviscous mantle case and p,T-dependent viscosity case
A.M. Bobrov, A.A.Baranov Schmidt Institute of Physics of the Earth, Russian Academy of Sciences, Moscow (bobrov@ifz.ru, baranov@ifz.ru)

In numerical two-dimensional experiments we investigate the spatial field of overlithostatic horizontal stresses in the mantle and in moving continent and its evolution. A continent moves self-consistently with changing mantle flows. Velocity of a continent in the process of movement varies in accordance with time-dependent forces which act from underlying viscous mantle as well as with mantle forces acting on the end faces of continent (i.e., the thickness of a continent is taken into account). This model is described in [Bobrov, Trubitsyn, 2008]. Continent viscosity is equal to 105 with respect to average viscosity of the mantle.

We consider two model laws for viscosity: isoviscous mantle case ν = const; and p,T-dependent viscosity case ν = 30 • exp [-9.2T + 2.3(1-z)] (z – depth from the surface). For these two models we analyze how a form of viscosity law can change the horizontal stress fields in the mantle and continent. We research what model law gives the results more close to actual data.

In our statement, both cases should give approximately equal Nusselt numbers (i.e., should have the same heat transfer). For this reason, the computational Rayleigh numbers were different (for exponential viscosity case - Ra = 107; for isoviscous model, where the heat transfer is more effective owing to the absence of high-viscosity layer at the surface - Ra = 106).

Numerical results
The horizontal stress field in moving continent greatly depends on variations of horizontal velocity in the underlying mantle, and also on continent position between the ascending and descending mantle streams. Sub-continental upwelling mantle flows have the extensive effect; sub-continental downwelling ones - the compressive effect. Mantle plumes near continent`s borders demonstrate compressive effect on continent, downwelling flows produce its extension. If the horizontal stresses are presented in dimensionless form then our two cases (constant and variable viscosity) show rather big but not principal differences (fig. 1a, b; note that the stress values in the case of variable viscosity are higher). Thus, the mantle model with constant viscosity can be regarded as qualitatively representative.

However, after transition to dimensional parameters it appears, that in isoviscous mantle model stresses values bigger in several times, than in the case of variable viscosity. In this case isoviscous mantle model leads to strongly overestimated stresses and is not representative in this aspect.

Mantle model with variable viscosity has typical horizontal stress values in the major portion of mantle – (2 ÷ 6) MPa; in continent at different stages of its movement – (2 ÷ 15) MPa. (Fig. 2 – three successive stages of the movement of a continent; dimensionless time - 2.2•10-4, 3.4•10-4, 4.6•10-4).

References
Bobrov A.M., Trubitsyn A.P. Numerical model of the supercontinental cycle stages: integral transfer of the oceanic crust material and mantle viscous shear stresses // Stud. Geophys. Geod. 2008. V. 52. P. 87-100.

Figure 1a, b. Dimensionless horizontal stress fields σ _xx (x,z) for isoviscous model and p,T-dependent viscosity model. Light grey colors correspond to positive values, dark grey – to negative values. Stresses σ _xx (x,z) are determined by the relation σ _xx (x,z)= p(x,z) – 2 ∂y_x(x,z) /∂x, i.e., compressive stresses are considered to be positive. White isolines show the temperature field.

Fig 2 a,b,c. Dimensionless horizontal stress fields for p,T-dependent viscosity (dimensionless time- 2•10^-4, 3•10^-4, 4•10^-4)

Top

Continental Tectonics, Rheological Stratification and Seismicity: Insights from geodynamic numerical modeling.
E. B. Burov ISTEP, University of Pierre et Marie Curie (University of Paris 6), 4 Place Jussieu, 75252 Paris Cedex 05, France

Depending on the loading and thermal conditions and time scale, the lithosphere exhibits elastic, brittle-plastic or viscous-ductile properties. As suggested by rock mechanics experiments, a large part of the long-term lithospheric strength is supported in the ductile regime. Unfortunately, these data cannot be reliably interpolated to geological time and spatial scales (strain rates ~10e-17 – 10e-13 1/s) without additional parameterization. An adequate parameterization has to be based on “real time” observations of large-scale deformation. For the oceanic lithosphere, the Goetze and Evan’s brittle-elastic-ductile yield strength envelopes derived from data of experimental rock mechanics were successfully validated by a number of geodynamic scale observations such as the observations of plate flexure and the associated Te estimates. For continents, the uncertainties of flexural models and of other data sources are stronger due to the complex stratified rheological structure and history of continental plates. For example, in one continental rheology model, dubbed “jelly sandwich”, the strength mainly resides in the crust and mantle, while in another, dubbed “crème-brûlée”, the mantle is weak and the strength is limited to the upper crust. These models have arisen because of conflicting results from earthquake, elastic thickness (Te) and rheology data. We address these problems here by reviewing rock mechanics data and by examining the plausibility of each rheological model from general physical considerations. We next review the elastic thickness (Te) estimates and their relationship to the seismogenic layer thickness (Ts). We then explore, by analytical and numerical thermo-mechanical modeling, the implications of a weak and strong mantle for tectonic structural styles. We show that, irrespective of the actual crustal strength, the “crémé-brûlée” model is unable to explain either the persistence of mountain ranges for long periods of time or the integrity of the down-going slab in collisional systems. We conclude that, despite its limitations, the “jelly sandwich” is a useful first-order model with which to parameterize the long-term strength of the continental lithosphere. We show that dry olivine rheology laws best represent the long term behavior of the mantle lithosphere (Figure1). As to the crustal rheology, analysis of Te derived from robust forward flexural models, and from strength estimates obtained from numerical thermo-mechanical models of tectonic deformation, suggests that most continental plates need a strong mantle and a weak lower or intermediate crust (Figure 1). The models based on the “jelly sandwich” rheology predict distribution of seismogenic zones that comply with the observed values of Ts, whereas those based on the “crème brulée” rheology tend to under-estimate Ts by a factor of 2-3. We suggest that “jelly-sandwich” rheology is more compatible with the observed Ts data and that the aseismic behavior of the mantle is related to increasing confining pressure (brittle strength), or to a change in the brittle law with depth, rather then to any “ductile weakness”. Consequently, seismicity is a sign of weakness, not strength. We conclude that the Ts has no direct relation to long term rheology but is indicative of the stress/strain state in the crust. Finally, we derive a classification of possible large-scale deformation styles (Figure 1 – example for convergent processes) as a function of thermo-rheological profile of the continental lithosphere.

Figure 1. Summary of the results of numerical thermo-mechanical visco-elasto-plastic models, showing the predicted continental collision styles as function of the tested thermo-rheological profile (yield-stress envelopes, YSE, in inserts below each model). Te is the equivalent elastic thickness of the lithosphere computed directly for each YSE. Tm is temperature at Moho depth. Top right insert shows an example of crustal-scale structures (faults, fold-and-belt structures, surface topography) predicted by the model for the YSE on the left. These predictions can be directly matched with the observations allowing for structural validation of the assumed crustal and mantle rheology.

Top

Slab Geometry and the Compatibility Condition for the Kinematics of Subduction
Ling-Yun Chiao ^1,2, Chen-Chien Peng²
¹ Institute of Oceanography, National Taiwan University
² Department of Geosciences, National Taiwan University

The kinematical field associated with the lithosphere subduction is essential for the understanding of arc volcanism, seismogenesis and paleo-plate reconstruction. There is however, no straightforward way to specify it because the geometry of the subducted slab is neither Euclidean nor is it analogous to the Earth’s surface with a nearly uniform Gaussian curvature. We show that the particular slab geometry of the Northwest Pacific subduction zone being a natural consequence of avoiding severe membrane deformation rates within the slab. We envision that the principle of the minimum deformation rates that defines a flow field on the non-Euclidean surface as a generalization of the classical surface plate kinematics. We also derive the general compatibility condition that characterizes the fundamental interplays among the local geometry (the Gaussian curvature), the kinematical field and the associated deformation rates on a general non-Euclidean surface,

To our knowledge, this general compatibility equation, for kinematics defined for the curvilinear coordinate specified on a general non-Euclidean surface, has not been discussed in the past. It specifies clearly that the variation of the Gaussian curvature manifested by the subduction flow field constrains the compatibility of the deformation rate of subduction flow.

Figure 1. (a) Slab geometry portrayed by 100km depth contours. (b) In-plane strain-rate tensor field described by the compressional axes (bold bars) and tensional axes (thin line segments) and (c) the subduction flow field (paths originated from the trench and subduct down-dip) overlaid upon the magnitude variation of the effective strain-rates (the gray scale is determined by the logarithms of the effective strain-rate, i.e., gray scale -14 to -15 is for effective strain rate with magnitude between 10-14 and 10-15 per second). The upper panels (a, b, c) are for the synthetic slab geometry that is constructed by first assigning subducting dip angles on cross sections marked by dash lines in the Mariana trench and the northern edge in Kuril trench in (a) and then linearly interpolating the geometry in between to examine the impacts without the arch structure in the real observation. The flow field is constructed by simply rotating the surface Eulerian kinematics onto the local slab surface (notice that it converges under Japan Sea, causing high in-plane strain-rates). The lower panels (d, e, f) are for presentations of the case that we held the prescribed dips only along selected cross sections, dash lines on (a), and invert for the slab geometry and the appropriate flow field on that surface that minimizes the integrated in-plane deformation rates. Notice how the resulting geometry mimics the actual slab revealed by the Wadati-Benioff seismicity, and how the flow field has been adjusted to avoid unreasonable in-plane deformation rates.

Top

General parameterization and resolution of geophysical inverse problems
Ling-Yun Chiao^1,2,*, Yuan-Cheng Gung², Ya-Ting Hsu1 Yu-Hsuan Chang², Hsin-Ying Yang²
¹Institute of Oceanography, National Taiwan University, Taipei, Taiwan, R.O.C.
²Department of Geosciences, National Taiwan University, Taipei, Taiwan, R.O.C.

It is becoming a progressing consensus in the practice of resolving geophysical inverse problems that to seek for a family of related inverse models rather than just a single solution would help significantly to shed lights on the essence of Earth structures that a particular data set is capable of robustly resolving. Conventionally, after an a priori finite parameterization by expanding in terms of an arbitrary choice of truncated basis functions, the family of solutions is constructed along a certain tradeoff curve, for example, between data variance reduction versus model roughness or model norm; or between quantitative measures of model resolution versus model variance. It is seldom explored that cross assessment of the consistency of inverse models invoking distinctive parameterizations and comparative appraisal across different regularization schemes are also essential to differentiate whether specific features in the obtained solution(s) might actually be artifacts arisen from the particular parameterization and the inverse algorithms invoked. Whereas most parameterizations are based on ad hoc choices of truncated bases functions, it is important to emphasize the role of the resolution kernel that relates the model estimates to the true Earth’s structure, and that it is always much more insightful than ambiguous impressions obtained from the usually improvised checkerboard test. Furthermore, it is worth mentioning that there are corresponding resolution kernels for every distinctive parameterization rather than just the usually much accentuated spatial resolution. Here we show that cross domain transformations of the inverse model itself, the corresponding discrete Frechet kernel, the resolution as well as the model covariance matrices can all be achieved directly by straightforward algebraic transformations rather than tediously re-implementing the forward data rule for every distinctive parameterization from the scratch. By taking advantage of these transformation rules, not only does the cross domain appraisal is made easy, it is also possible to formulate the forward data rule in one basis such as a global function set and then to examine the inverse models in terms of the other basis such as a local set and vice versa. This offers great flexibility to take advantage of the intrinsic merits of each of every distinctive basis functions. In short, straightforward avenue for inspecting a complete solution, including the inverse model and the corresponding model resolution and covariance matrices, in terms of different parameterization basis, different weighting and/or different regularization schemes becomes immediately available. That is, the sometimes appealing desire, while interpreting the implications of a certain inverse model, to reiterate an existing suit of solutions by evaluating the effect of reweighting or reparameterization can be done through simple algebraic transformation rather than computing the usually sizeable matrix inversion over and over again. The transformation rule is summarized as,

Figure 1. The results of the multiscale based inversion (upper row) as compared with the initial global model SAW642AN_iso (Panning and Romanowicz, 2006). An additional 25% variance reduction is obtained for the optimal solution chosen along the tradeoff curve.

Top

Lattice thermal conductivity of MgO at conditions of Earth’s interior
Jianjun Dong^a and Xiaoli Tang^a,b
a. Physics Department, 206 Allison Laboratory, Auburn University, Auburn, AL 36849
b. Earth and Space Science Department, University of California, Los Angeles, 595 Charles Young Drive East, Box 951567, Los Angeles, CA 90095.

We applied a recently developed first-principles computational method to predict the lattice thermal conductivity of (κ_MgO) over a wide range of temperature-pressure conditions within the framework of Peierls-Boltzmann transport theory. A density and temperature dependent function of the thermal conductivity for MgO (κ_MgO(ρ,T)) is proposed for conditions from ambient to the Earth’s core-mantle-boundary (CMB). κ_MgO at lower mantle conditions was then calculated based on the proposed κ_MgO(ρ,T) model and the LDA calculated thermal equation of state . Our calculated κ_MgO and its pressure derivative d lnκ/dP at ambient condition are 66 W/K/m and 3.9% GPa^-1 respectively, which are in good agreement with experiments. At the mantle side of the CMB (roughly T=3000K, P=135GPa), we predict κ_MgO and d lnκ_MgO/dP to be 43.4 W/K/m and 0.36% GPa-1 respectively, suggesting MgO contributes about 8.7 W/K/m to the total thermal conductivity at the top of the CMB. Furthermore, our calculation shows that κ_MgO increases with depth, and falls in the range between

and

based on the hot and cold geotherms constrained by experiments.

Top

Wave-packets and their applications in seismology
Duchkov A.A. ¹, Andersson F.², De Hoop M.V.³ 1 – Purdue University, West Lafayette, and IPGG SR RAS, Novosibirsk, aduchkov@purdue.edu
2 – Purdue University, West Lafayette, and Lund University, Lund, aduchkov@purdue.edu
3 – Purdue University, West Lafayette, mdehoop@purdue.edu

Introduction
In this talk we introduce wave packets, describe their properties and mention briefly their applications in seismology. We discuss a multi-scale geometric representation of (seismic) waves via decomposition into wave packets. Wave packets can be thought of as certain localized “fat” plane waves. There are different ways to introduce a discrete frame wave packets for decomposing data. One such discretization leads to the introduction of digital curvelets [2]. Here, we follow a different discretization while constructing discrete almost symmetric 3D wave packets (this is a generalization of a 2D frame introduced in [1]).

We get a discrete transform for decomposition-recomposition of a 3D data set into wave packets. Wave-packet representation of seismic data appears to be sparse, i.e. comparatively small amount of wave packets provides a good approximation of the data. Our transform is unitary. Another relevant aspect of our discretization is the appearance of parameters that control the form of wave packets used in the transform. This form can be adapted to the physical problem at hand. We consider applications in exploration seismology and global seismology. Possible applications of the multi-scale, discrete transform developed here, include higher-dimensional data regularization (including l1 optimization), directional regularity analysis, feature extraction and reflection tomography.

Properties of wave packets
Wave packets being oscillatory in one direction and smoothly varying in other directions can be easily associated with local plane waves, i.e. pieces of a wavefront. Every wave packet is characterized by its central pint and orientation (direction in which it is oscillatory). Thus it appears natural that wave-packet decomposition with subsequent thresholding of coefficients extracts a comparatively small number of wave packets. These packets are clustered with central points being close to traveltime surfaces and orientations being close to normal directions of these traveltime surfaces. Thus one can perform approximate detection of seismic events by simply transforming the data. Another property of wave packets is that they are, in some sense, almost eigen-functions for solving initial value problems for a wave equation or other hyperbolic evolution equations. It means that initial data can be decomposed into wave packets and these packets can be simply translated along a Hamilton flow associated with an evolution equation. As a result we get a good approximation to a solution of the initial problem that can be further corrected by solving a corresponding Volterra integral equation.

Wave-ray duality
Wave packets provide a connection between ray-based kinematic methods and wave-equation based methods designed for band-limited signals. The connection is based on the wave-ray duality of wave packets. On the one hand wave-packet frame provides numerically exact decomposition-recomposition of band-limited seismic data. On the other hand after decomposition each coefficient can be identified with corresponding wave packet which has specific geometric characteristics: position of its central point and orientation.

Reflection tomography
We use wave packets for the purpose of reflection tomography, i.e. finding a smooth background velocity model from reflected waves present in data. Compared to a transmission tomographic inversion, in which case source and receiver locations are usually known, in reflection tomography a reflector position is another unknown. Thus it is not a surprise that redundancy in seismic data (data set dimension should be greater then subsurface domain dimension) is needed for successful reflection tomographic inversion. This redundancy is usually exploited in form of optimization function that measures a semblance criterion along a specific parameter an specially constructed gathers: angle gathers or common-image gathers (cf. [4]).

Wave packets can be used to represent wave-equation derived common-image gathers (cf. [4]) as well as provide geometrical information that can be further utilized in “ray”-geometrical analysis, which is exploited in the context of reflection tomography. We propose a new semblance criterion for the common-image gathers based on the orientation of wave packets (as an alternative to the residual moveout). Then we apply linear perturbation theory to the geometry of the angle transform in order to derive formulas for the background velocity update.

Conclusions
Discrete almost symmetric wave packets are well suited to develop a wave-ray duality in computational propagation and scattering problems. This duality is apparent in the clustering of wave-packet coefficients of wavefields and corresponding compression following wave fronts in seismic data sets. Also wave packets allow “ray”-geometrical formulation of a reflection tomography inverse problem for finding smooth background velocity model from reflected waves in seismic data.

References
1. Andersson F., De Hoop M.V., Smith H., Uhlmann G. A multi-scale approach to hyperbolic evolution equations with limited smoothness // Comm. Partial Differential Equations. 2008. Vol. 33, pp. 988–1017.
2. Candes E., Demanet L., Donoho D., Ying L. Fast discrete curvelet transforms // SIAM Multiscale Model. Simul. 2006. Vol. 5-3, pp. 861–899.
3. Stolk C.C. and De Hoop M.V. Seismic inverse scattering in the downward continuation approach // Wave Motion. 2006. Vol. 43, pp. 579 – 598.
4. Symes W.W., Carazzone J. Velocity inversion by differential semblance optimization // Geophysics. 1991. Vol. 56, pp. 654 – 663.

Top

Constraints on rheology of plastoshere from space geodetic observations of deformation following large earthquakes
Yuri Fialko Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography University of California, San Diego La Jolla, CA 92093-0225

Large shallow earthquakes generate sudden stress perturbations in the ambient lithosphere that have magnitude of order of megapascals, and spatial wavelength of order of 104-105 meters, i.e., sufficiently large to induce a ductile response of the lower crust and the uppermost mantle. Much improved global imaging of surface deformation due to the existing Interferometric Synthetic Aperture Radar (InSAR) satellites allow us to measure the delayed postseismic response of the crust and the upper mantle to the coseismically induced stress changes in areas of recent large earthquakes. Measurements of postseismic transients bear on a long-standing debate on the effective rheology of the lithosphere below the brittle-ductile transition. Commonly considered models of postseismic transients include enhanced creep on a seismic rupture or its extension below the brittle-ductile transition (the so-called afterslip), poro-elastic rebound of the fluid-saturated crust, and visco-elastic relaxation of the lower crust or upper mantle. A robust discrimination (or some assessment of relative importance) of these mechanisms is often hindered by similarities in the predicted patterns of surface displacements. Such discrimination is nevertheless critical for our understanding of post-earthquake stress changes, delayed triggering of seismicity, and long-term strength of the ductile part of the continental lithosphere. I will present space geodetic (InSAR and GPS) observations of postseismic transients from several large earthquakes in tectonically diverse regions, including the 1992 M7.3 Landers and 1999 M7.1 Hector Mine events in the Eastern California Shear Zone, the 2003 M7.2 Altai event at the northern boundary of the Tibetan plateau, the 2006 M7.0 Mozambique earthquake at the southern tip of the East African rift zone, and the 2008 M7.9 Sichuan earthquake at the transpressional eastern boundary of the Tibetan plateau. I will argue that the data are inconsistent with models of low-viscosity and low-strength lower crust (e.g., Figure 1). The observed deformation signals (or in some instances, the lack of observable signals) suggest that the effective strength of the lower crust is higher than that of the upper mantle. The essentially non-exponential temporal decay of the observed transients cannot be explained in the framework of simple linear viscoelastic models; possible alternatives include afterslip (consistent with laboratory rate-and-state phenomenology), power-law creep and bi- (or multi-) viscous rheology. Spatial patterns of postseismic displacements provide further constraints on admissible constitutive relations. While it is recognized that deep afterslip and viscoelastic relaxation are kinematically indistinguishable in case of two-dimensional deformation (e.g., under anti-plane strain conditions relevant to very long strike-slip ruptures), recent work has shown that the surface displacement patterns predicted in case of 3-D loading (e.g., dipping faults, rupture sizes of order of thickness of the brittle layer) may allow one to discriminate between possible mechanisms, provided that both horizontal and vertical components of the displacement field are available. Also, the focus has shifted from kinematic inversions for an arbitrary distribution of fault slip to physically-based models of afterslip in which the latter is driven by coseismic stress changes on a fault plane (thereby preventing excessive slip in areas that didn't experience much loading by an earthquake). The degree of strain localization below major seismogenic faults (the so-called “fault roots”) critically affects the spatiotemporal patterns of postseismic transients.

Existing models are often limited in their ability to reproduce the relevant physics involved in postseismic relaxation. I will discuss available numerical approaches that allow us to simulate all three phenomena believed to be important in postseismic response: stress-driven afterslip, non-linear viscoelastic relaxation in the bulk of the ductile substrate, and time-dependent poroelasticity.

Figure 1. Post-seismic interferograms and predicted viscoelastic relaxation from the 2003 Mw7.2 Altai event. (a) Radar interferogram spanning 2003/10/13-2006/08/28, ENVISAT track 391. (b) Interferogram spanning 2003/11/01-2006/09/16, ENVISAT track 162. (c) Predicted viscoelastic relaxation in the lower crust (assumed Maxwell relaxation time of 0.7 years). (d) Predicted viscoelastic relaxation in the upper mantle (same Maxwell time as in (c)). Note a disagreement between observations (a-b) and model predictions assuming a weak lower crust (c).

Top

Use of geoinformation models for the analysis of stress-strain state of the earth's crust of the Caspian region
I.A. Garagash, A.V. Dubovskaya Institute of the Earth Physics, Moscow, Russia

The stress state of the earth's crust is far from lithostatic. This is caused primarily by the movements of lithosphere plates, which generate significant tangential stresses and high level of seismic activity close to the plate borders. Consequently, knowledge of stress distributions is necessary for seismic zoning and determining the location of possible strong earthquakes. On the other hand, non-uniformity of stress distribution in the earth's crust controls it fissuring and fluid flow, affecting formation of hydrocarbon deposits. It is therefore apparent, that robust quantitative estimations of stress distributions in different areas of the earth’s crust are needed.

Figure 1. Model of the Caspian region.

Modern computing as well as the level of geological and geophysical information available about the earth's crust and lithosphere allow creation of 3D geomechanical models for many places on the Earth. Such a 3D model, including the basic structural borders, is created here for the Caspian region, using the geophysical data . The external loads generating the stress state are the body weight, temperature and the horizontal tectonic movements. The latter are entered into calculation on the basis of global velocity distribution model of plate movement from the GPS data.

The model of the earth's crust of the Caspian region is created for the area from 37.0°N – 47.5°N to 46.0°E – 56.0°E. The model contains the mountain relief, bathymetry, and characteristic internal borders such as the roof of the consolidated crust, Moho surface and upper asthenosphere boundary (Fig. 1). The characteristic size of the studied area is 870kmX1200kmX180km. The asthenosphere surface for the Caspian region has been received on the basis of the heat flux data in region, taking into account the internal structure of crust and thermal properties of rocks.

It is assumed that rocks behave as an elasto-plastic material with Drucker-Prager yield condition. Asthenosphere is modeled by the viscoelastic Maxwell material. Conditions of unrestricted vertical sliding are prescribed along the vertical boundaries of the model.

The calculation is done in two stages. First the initial stress state of the model under action of gravitational forces and temperature distribution is calculated. Then the loading by the tangential tectonic forces generated by the tectonic plate movements is applied. For the latter purpose we use the standard model of velocity distribution plate movements NNR-NUVEL.

As a result of the calculation, the data about stress state invariants distribution in the region are obtained. These parameters, along with the computed distribution of the accumulated elastic energy, allow for the seismic zoning of the region. This work was funded by grant RFBR №08-05-01027.

Top

DEPTH ESTIMATION FOR THE BASEMENT SURFASE USING MODIFIED SPECTRAL ANALYSIS TECHNIQUE
G.S.Hassan Minia University, Faculty of Engineering, Petroleum Engineering Dept. E-mail: g seliem @ yahoo. com

Basement surface is one of the most important targets that can be determined from the gravity data in order to interpret adequately the structure geology of the area. Moreover, the depth to the top of the basement complex can give the thickness variations of the sedimentary section overlying the basement surface.

Most of researchers used modulus logarithm of gravity and magnetic anomaly spectra to determine the upper (Z1) and lower (Z2) surfaces of causative masses, but without corrections for these depths.

The developed formulae for spectral analysis technique which we suggested in many published papers represent an important tool to determine the average depth values to buried bodies precisely. The depth of the causative geological bodies is the most important parameters that can be determined form gravity and magnetic anomalies.

Appling the present modified technique of spectral analysis to gravity and magnetic anomalies of calculating the depth to the basement surface configuration and comparing these results with the other geophysical methods showed accurate results.

Top

Crustal deformation deduced from micro-gravity and geodetic data in the High Dam area, Aswan (Egypt).
G.S. Hassan Egypt, Minia University, Faculty of Engineering, Petroleum Engineering Dept. E-mail: gseliem @ yahoo.com

After the 1981, November 14 th earthquake, the Aswan High Dam area in Egypt has been subject to various geophysical and geodetic studies, due to its importance for some vital projects. Gravity studies were performed in the form of detailed measurements, in order to determine the basement relief in this area.

In addition, the micro-gravity measurements were accomplished by the GPS and precise leveling measurements, in order to determine the accurate height of the observation points. The resulted gravity data were adjusted and processed and the Bouguer anomaly map of the area was obtained. The derived depths to the basement are ranging between few meters and 200 meters. Furthermore, good information about the distribution of the fault system in the study area was concluded.

The type of the prevailing faults in the area is normal faulting, while the main trends of faults are the N-S, E-W and NE –SW directions. The accuracy of the GPS and leveling measurements is of 23 cm. The low level of earthquake occurrences and the low strain rates indicate that, the rate

Top

Scale aspects of heat transport in the diamond anvil cell, in spectroscopic modeling, and in Earth’s mantle
Anne M. Hofmeister, Department of Earth and Planetary Sciences, Washington University, St. Louis, Missouri 63130-4899, U.S.A. [hofmeist@levee.wustl.edu; (314) 935-7440 fax: (314) 935-7361]

Examining length and time scales reveals severe problems in experimental and spectroscopic determinations of thermal conductivity (k) involving the diamond anvil cell. Laboratory conditions are optically thick (diffusive) for low frequency phonons but optically thin (ballistic) for high frequency photons. Because speeds differ by a factor of 105, photons carry all the heat ballistically when these processes operate in parallel, as occurs in highpressure,transient heating experiments recently published (e.g. Beck et al. J.Appl Phys. 2007; Goncharov et al. PEPI 2009) Not recognizing this effect lead to inappropriate processing of portions of the cooling curves as vibrational transfer when instead these depict ballistic radiative transfer from the embedded metal foil to the surroundings, thus providing bogus phononic k-values for MgO and other solids at lower mantle pressures. In spectroscopic experiments, scattering from imperfections and refraction effects involving the diamonds have been misinterpreted as sample absorption (e.g. Goncharov et al. Science 2006) and in addition Keppler et al. (Science 2007) used samples too thin to establish Fe2+ band strengths, leading to erroneous photonic k-values for lower mantle phases. Improved estimates of lower mantle k are provided by re-analyzing available data.

Scaling aspects furthermore raise questions as to how radiative transfer in the deepest mantle should be handled, and offer perspective on the importance of layering to cooling of the Earth.

Top

Three-dimensional numerical modeling of crustal extension
R.S. Huismans, V. Allken, and C. Thieulot Dep. Earth Science, University of Bergen, Norway (Ritske.Huismans@geo.uib.no)

Here we focus on understanding the evolution and structural style of crustal extension in 3D using state of the art computational modeling techniques. To date few 3D models exist that follow the evolution of tectonic processes into large deformation modes with sufficient resolution to resolve individual faults and shear zones. We use an Arbitrary Lagrangian Eulerian (ALE) thermo-mechanically coupled fully parallel Finite Element code which solves for visco-plastic flows in 3D. Plastic materials weaken with accumulating strain. To localize deformation, a weak seed region is introduced at the base of a one layer model extended by velocity boundary conditions. Controls on the geometry and spacing of three-dimensional frictional-plastic shear zones in simple one and two-layer models are investigated. We specifically focus on factors controlling rift propagation and rift obliquity for varying weak seed extent and obliquity with respect to far-field extensional boundary conditions.

Top

Data assimilation in problems of mantle dynamics
Alik Ismail-Zadeh1^2,3, Gerald Schubert ⁴, Igor Tsepelev ⁵, Alexander Korotkii ⁵
1 Geophysikalishes Institut, Universität Karlsruhe, Karlsruhe, GERMANY; e-mail: alik.ismail-zadeh@gpi.uka.de
2 MITPAN, Russian Academy of Sciences, Moscow, RUSSIA; e-mail: aismail@mitp.ru
3 Institut de Physique du Globe de Paris, FRANCE; e-mail: aiz@ipgp.jussieu.fr
4 Department of Earth and Space Sciences, University of California, Los Angeles, USA; e-mail: schubert@ucla.edu
5 Institute of Mathematics and Mechanics, Russian Academy of Sciences, Yekaterinburg, RUSSIA; e-mail: tsepelev@imm.uran.ru; korotkii@imm.uran.ru

To restore thermal structures in the mantle (e.g., ascending plumes and descending lithospheric plates) in the geological past, data assimilation can be used to constrain the initial conditions for the mantle temperature and velocity from their present observations. The initial conditions so obtained can then be used to run forward models of mantle dynamics to restore the evolution of mantle structures. Data assimilation can be defined as the incorporation of observations (in the present) and initial conditions (in the past) in an explicit dynamic model to provide time continuity and coupling among the physical fields (e.g., velocity, temperature). The basic principle of data assimilation is to consider the initial condition as a control variable and to optimize the initial condition in order to minimize the discrepancy between the observations and the solution of the model.

If heat diffusion is neglected, the present mantle temperature and flow can be assimilated using the backward advection (BAD) into the past. Two- and three-dimensional numerical approaches to the solution of the inverse problem of the Rayleigh-Taylor instability were developed for a dynamic restoration of diapiric structures to their earlier stages (e.g., [1-4]). The mantle flow was modeled backwards in time from present-day mantle density heterogeneities inferred from seismic observations (e.g., [5, 6]). Both direct (forward in time) and inverse (backward in time) problems of the heat (density) advection are well-posed. This is because the time-dependent advection equation has the same form of characteristics for the direct and inverse velocity field (the vector velocity reverses its direction, when time is reversed). Therefore, numerical algorithms used to solve the direct problem of the gravitational instability can also be used in studies of the time-reverse problems by replacing positive time-steps with negative ones.

The variational (VAR) (or also called adjoint) data assimilation has been pioneered by meteorologists and widely used in oceanography and in hydrological studies. The use of VAR data assimilation in models of geodynamics (to estimate mantle temperature and flow in the geological past) has been put forward by Bunge et al. [7] and Ismail-Zadeh et al. [8, 9] independently in 2003. The VAR approach by Ismail-Zadeh et al. [9] is computationally less expensive, because it does not involve the Stokes equation into the iterations between the direct and adjoint problems, and this approach admits the use of temperature-dependent viscosity. The VAR data assimilation algorithm was employed to restore numerically models of present prominent mantle plumes to their past stages [10] and to recover the structure of mantle plumes prominent in the past from that of present plumes weakened by thermal diffusion [11]. The VAR method was recently used to study dynamics models of thermal plumes and lithospheric plates in the mantle (e.g., [12-14]).

The use of the quasi-reversibility (QRV) technique implies the introduction into the backward heat equation of the additional term involving the product of a small regularization parameter and a higher order temperature derivative. The data assimilation in this case is based on a search of the best fit between the forecast model state and the observations by minimizing the regularization parameter. The modified QRV method was recently introduced in geodynamic modeling and employed to assimilate data in models of mantle dynamics [15, 16]. The advances in numerical modeling and in data assimilation attract an interest of the community dealing with dynamics of the mantle structures. The aim of this talk is to review the data assimilation methods used in geodynamic modelling and to present applications of these methods to geodynamical problems.

1. Ismail-Zadeh AT, Talbot CJ, Volozh YA (2001) Dynamic restoration of profiles across diapiric salt structures: numerical approach and its applications. Tectonophysics 337: 21-36.
2. Kaus BJP, Podladchikov YY (2001) Forward and reverse modeling of the three-dimensional viscous Rayleigh-Taylor instability. Geophys Res Lett 28: 1095-1098.
3. Korotkii AI, Tsepelev IA, Ismail-Zadeh AT, Naimark BM (2002). Three-dimensional backward modeling in problems of Rayleigh-Taylor instability. Proc Ural State Univ 22: 96-104.
4. Ismail-Zadeh AT, Tsepelev IA, Talbot CJ, Korotkii AI (2004) Three-dimensional forward and backward modelling of diapirism. Tectonophysics 387: 81-103.
5. Steinberger B, O’Connell RJ (1998) Advection of plumes in mantle flow: implications for hotspot motion, mantle viscosity and plume distribution. Geophys J Int 132: 412-434.
6. Conrad CP, Gurnis M (2003) Seismic tomography, surface uplift, and the breakup of Gondwanaland: Integrating mantle convection backwards in time. Geochem Geophys Geosystems 4: doi:10.1029/2001GC000299.
7. Bunge HP, Hagelberg CR, Travis BJ (2003) Mantle circulation models with variational data assimilation: Inferring past mantle flow and structure from plate motion histories and seismic tomography. Geophys J Int 152: 280-301.
8. Ismail-Zadeh AT, Korotkii AI, Tsepelev IA (2003) Numerical approach to solving problems of slow viscous flow backwards in time. In: Bathe KJ (ed) Computational Fluid and Solid Mechanics, Amsterdam: Elsevier Science, 938-941.
9. Ismail-Zadeh AT, Korotkii AI, Naimark BM, Tsepelev IA (2003) Three-dimensional numerical simulation of the inverse problem of thermal convection. Comput Math & Math Phys 43: 587-599.
10. Ismail-Zadeh A, Schubert G, Tsepelev I, Korotkii A (2004) Inverse problem of thermal convection: Numerical approach and application to mantle plume restoration. Phys Earth Planet Inter 145: 99-114.
11. Ismail-Zadeh A, Schubert G, Tsepelev I, Korotkii A (2006) Three-dimensional forward and backward numerical modeling of mantle plume evolution: Effects of thermal diffusion. J Geophys Res 111: doi:10.1029/2005JB003782.
12. Hier-Majumder CA, Belanger E, DeRosier S, Yuen DA, Vincent AP (2005) Data assimilation for plume models. Nonlinear Processes in Geophysics 12: 257-267.
13. Liu L, Gurnis M (2008) Simultaneous inversion of mantle properties and initial conditions using an adjoint of mantle convection. J Geophys Res 113, doi:10.1029/2008JB005594.
14. Liu L, Spasojevic S, Gurnis M (2008) Reconstructing Farallon plate subduction beneath North America back to the Late Cretaceous. Science 322: 934-938.
15. Ismail-Zadeh A, Korotkii A, Schubert G, Tsepelev I (2007) Quasi-reversibility method for data assimilation in models of mantle dynamics. Geophys J Int 170: 1381-1398.
16. Ismail-Zadeh A, Schubert G, Tsepelev I, Korotkii A (2008) Thermal evolution and geometry of the descending lithosphere beneath the SE-Carpathians: An insight from the past. Earth Planet Sci Lett 273: 68-79.

Top

Constraining the effects and amount of serpentinization in the oceanic lithosphere
Karthik Iyer and Lars Ruepke IFM-GEOMAR, Wischhofstrasse 1-3, 24148 Kiel, Germany. *Email: kiyer@ifm-geomar.de Tel: +49 431 600 2866 Fax: +49 431 600 2941

Seawater circulation through the oceanic lithosphere and associated hydration reactions play a key role in many geological processes and settings. Two of these settings that may significantly impact the amount of water assimilated into the mantle are mid-ocean ridges and subduction zones. A heat source in the form of a magmatic body at the ridge as well as numerous tectonic features such as faults drives hydrothermal circulation and brings fluids down to the mantle. At active margins and subduction zones, bending and faulting provide pathways through which fluids may reach the underlying mantle by seismic pumping. One key hydration reaction is the serpentinization of mantle rocks. The transformation of a ‘dry’ peridotite to a ‘wet’ serpentinite results in an uptake of ~13wt.% of water, a volume increase and density decrease of ~40%, and a strong decrease in mechanical strength. These drastic changes in rock properties illustrate the importance of quantifying the degree of serpentinization caused by deep seawater circulation. Furthermore, the amount of water stored in the mantle which is later recycled may also be quantified. For this purpose, we have developed a new hydrothermal convection model that also accounts for serpentinization. The key feedbacks of the reaction on fluid flow, i.e. variations in porosity/permeability due to volume changes, exothermic heat of reaction and reaction-induced fluid consumption, are all accounted for. We have applied this model to two test cases: hydrothermal convection at mid-ocean ridges and bend-faulting related hydration of subducting plates. For the mid-ocean ridge case, we have coupled a kinematic thermal ridge model to the hydrothermal convection model. The coupled model allows us to study the feedbacks between serpentinization and hydrothermal flow. We find that, on the one hand, hydrothermal convection is not confined to the crust but may extend well below the Moho. On the other hand, extensive serpentinization has the ability to close pore space and reduce permeability thereby hindering flow. This serpentinized layer, therefore, acts as a self-limiting boundary to hydrothermal flow. Additionally, the amount of serpentinization occurring at mid-ocean ridges is inversely related to the spreading rate giving rise to two distinct regimes for water content at slow- and fast-spreading ridges. At the other end of the lithospheric cycle, bend-faulting of subducting lithosphere, may provide the pathways for seawater to reach and react with the cold lithospheric mantle to make serpentine. To test this hypothesis, we explore the conditions under which seawater may reach the mantle and what the likely hydration pattern around fault zones are. Seismic velocities are related to the extent of serpentinization in the mantle and can therefore be compared to measurements taken at subduction zones (e.g. Middle-America trench).

Top

Influence of pressure on the deformation fabrics of olivine
H. Jung^1*, W. Mo1, H. W. Green², M. Park¹, J. Lee1, and S. Jung¹
1. School of Earth and Environmental Sciences, Seoul National University, Seoul 151-747, Korea (hjung@snu.ac.kr).
2. Institute of Geophysics and Planetary Physics and Department of Earth Sciences, University of California, Riverside, CA 92521, U. S. A.

Seismic anisotropy (SA) in Earth’s upper mantle is widely observed and considered to be caused by deformation-induced lattice-preferred orientation (LPO) of olivine. SA within oceanic lithosphere is attributed to LPO of dry olivine induced by large strains during corner flow at ocean ridges. More recent work shows that flow under high H2O fugacity induces a change in olivine LPO that explains the SA in the mantle wedge above subducting lithospheric slabs. Whether changes in olivine LPO are unique to fH2O effects has become controversial and is critical to resolve. Here, we report low-stress, high strain, experiments on dry harzburgite (96% olivine) at T = 1300°C and P = 2.5-3.6GPa. We show that at ~3GPa, pressure induces the same profound transition in olivine LPO that is produced by high fH2O at 1 – 2 GPa. One important consequence for global tectonics is that trench-parallel SA of the fast S-wave beneath subducting slabs probably reflects entrainment of asthenospheric mantle beneath the slab in the direction of subduction rather than trench-parallel flow as currently interpreted. The variety of olivine LPOs in both experiments and natural rocks suggest that, in addition to the pressure-induced change in olivine slip systems implied here, there are likely further changes in slip systems at higher pressure and temperature. Recent studies on the petrofabrics of olivine in mantle xenoliths will also be discussed.

Top

Some issues on rheological properties of minerals pertinent to global geodynamics
Shun-ichiro Karato Yale University, Department of Geology and Geophysics

Rheological properties of minerals control the dynamics of terrestrial planets. However, experimental studies on rheological properties are challenging not only because well-controlled deformation experiments under deep planetary conditions are difficult but also because of the large number of factors that could affect the rheological properties. In this talk, I will discuss a few topics related to the rheological properties of minerals that have a close link to the global geodynamics:

(1) The strength of the lithosphere: why plate tectonics on Earth? Previous models for lithosphere strength were based on the properties of olivine. Our recent experimental studies show that orthopyroxene has much smaller creep strength under the lithospheric conditions. I will discuss the possible implications of these observations for the operation of plate tectonics on Earth but not on Venus.

(2) The energy dissipation in convecting mantle: some notes on the scaling law for thermal evolution of terrestrial planets Thermal evolution of terrestrial planets have been studied using a “parameterized convection” model where some relationship between the average temperature and energy transport by convection are used. In a conventional model, the rate of energy transport increases with temperature, but Conrad-Hager (1999, 2000) (also Korenaga, 2003) suggested that when energy dissipation of plate bending plays an important role, then such a relation might need to be modified. I will discuss the validity of such a modification and present an alternative model for energy dissipation based on the latest knowledge of rheological properties of minerals under deep mantle conditions.

Top

Effective Viscosity of Suspensions of Swimming Bacteria
Dmitry A. Karpeev Argonne National Laboratory and Computation Institute University of Chicago, Chicago, IL karpeev@mcs.anl.gov

We present a model for dilute suspensions of swimming bacteria in a Stokesian uid and use it to calculate the eective viscosity of the suspen- sion. A bacterium rotates its agella to propel itself forward or to reorient itself randomly in a process called tumbling. Due to a bacterium's asymmetric shape, interactions with a prescribed generic background flow, such as a purely straining or a planar shear flow, cause the bacterium to preferentially align in certain directions. In these aligned configurations, self-propulsion produces a reduction in the effective viscosity. For suficiently weak background flows, the effect of self-propulsion on the eective viscosity dominates all other contributions, leading to an eective viscosity of the suspension that is lower than the viscosity of the ambient fluid and, potentially, even a locally negative effective viscosity. This is in agreement with recent experiments on suspensions of Bacillus subtilis. We present a number of asymptotic results for the dilute limit, in which the eective viscosity can be explicitly quantified in terms of the geometry of the bacteria, the parameters of self-propulsion and the characteristics of the background flow. This is joint work with Leonid Berlyand, Brian Haines, Vitaliy Gyrya (Penn State) and Igor Aronson (Argonne/U.Chicago).

Top

Rheology, deformation, shear localization and strength of the lithosphere.
B. Kaus¹, S. Schmalholz², F. Crameri¹, J.-P. Burg²
(1) Institute of Geophysics, ETH-Zurich, Sonnegstrasse, 5, 8092 Zurich, Switzerland.
(2) Institute of Geology, ETH-Zurich, Sonnegstrasse, 5, 8092 Zurich, Switzerland.

The understanding of thermo-mechanical processes forming mountain ranges and plateaus, such as the Tibetan Plateau, is still incomplete. Therefore, interpretations of the same geophysical observations with respect to lithospheric processes in general and lithospheric rheology in particular vary signiﬁcantly. This is, for example, the case for interpretations of a weak lithospheric mantle rheology based on the relative scarcity of earthquakes in the sub-continental mantle, where lithospheric strength proﬁles generally predict extremely high differential stresses. The weak mantle assumption, however, implicitly assumed that strength envelopes are valid in actively deforming lithospheric regions.

We test this assumption on two end-member lithospheres with (1) a weak lower crust and strong mantle and (2) a strong lower crust and weak mantle. For this purpose, we compare 1D models with 2D visco-elasto-plastic numerical models of continental shortening. Both 2D models show that strongly heterogeneous deformation typically follows initially homogeneous deformation. Lithospheric-scale buckle folds and shear zones result in strain rate variations of up to three orders of magnitude. Differential stresses in the upper crust are close to the yield strength, as predicted by 1D models. Stresses in deeper lithospheric regions, however, are signiﬁcantly smaller than in 1D models, especially in actively deforming regions [1].

Systematic numerical simulations as a function of temperature and deformation rate reveal that 1D models are reliable in hot and/or slowly deforming lithospheres only. The relative scarcity of earthquakes at mantle levels should thus be interpreted as an intrinsic consequence of strong lithospheric deformation rather than as evidence for a weak upper mantle rheology.

The simulations also highlight the importance of lithospheric-scale buckling and shear localization. To further understand the mechanisms that are responsible for lithospheric-scale shear-localization, we have performed a large number of 2D visco-elasto-plastic simulations. Results indicate that shear-heating, in combination with large differential stresses and low-temperature plasticity is an important factor to initiate localization. We are currently in the process of comparing the 2D results with the results of 1D simulations and previously derived scaling laws for the onset of shear-heating induced localization in viscoelastoplastic rocks [2], with the aim to derive and further unravel the key parameters that may control lithospheric failure.

References.
[1] Schmalholz S.M., Kaus B.J.P., Burg J.-P. (in press). Stress - strength relationship in the lithosphere during continental collision. Geology.
[2] Kaus B.J.P., Podladchikov Y.Y. (2006) Initiation of localized shear zones in viscoelastoplastic rocks. Journal of Geophysical Research. Vol 111, B04412, doi: 10.1029/2005JB003652.

Top

Classical Density Functional Theory for Hard-Sphere Liquids
Matthew G. Knepley Computation Institute, University of Chicago, Chicago, IL knepley@ci.uchicago.edu May 27, 2009

Classical Density Functional Theory (DFT) [4] provides a theoretical framework for calculation of a electrochemical potential for charged hard spheres in a dielectric liquid. This framework has recently been used to simulate the permeation of ions through protein channels [2, 3]. However, DFT can address the more general problem of hard spheres in suspension, accomodating several species of varying radii. We may also incorporate long range correlations using a method similar to the treatment of electrostatics in [2], in which the chemical potential of the long range interaction is expressed as an expansion around a local equilibrium solution, which in the case of electrostatics is the Mean Spherical Approximation (MSA). In its current formulation, DFT can describe a steady flow when used in a Nernst-Planck formulation, which may be appropriate for slow geophysical flows. For unsteady flows, it may be possible to employ a recent dynamic scheme [1]. A generic DFT package should therefore allow rapid exploration of novel rheologies arising from geophysical solid-liquid mixtures, such as slurries, melts, and perhaps some magmas. The computational structure of our current implementation, based upon PETSc, will be discussed.

Top

Мегациклический режим термохимической конвекции земной мантии
Котелкин В.Д. ¹, Лобковский Л.И. ²
¹ - МГУим. М.В.Ломоносова, Москва, kotelkin55@mail.ru
² - ИОРАН им. П.П.Ширшова, Москва, llobkovsky@ocean.ru

Исследования мантийной конвекции с фазовыми переходами неоднократно проводились [1-3] и проводятся [4] на основе термической модели. Установлено, что эндотермический фазовый переход тормозит термическую конвекцию и ведет к перемежающимся режимам. Альтернативная отечественная химико-концентрационная схема конвекции [5, 6] долгое время не находила поддержки и развития. Однако на фактическом материале было показано [7], что мантийная динамика реализуется двумя способами: непрерывно действующей тектоникой плит и импульсно действующей тектоникой плюмов. Термическая конвекция в состоянии объяснить только тектонику плит, а внутриплитные объекты – горячие точки, траппы и океанические поднятия – не находят объяснения в её рамках. Мы, вот уже более 10 лет, проводим численные исследования, объединяя оба подхода в рамках термохимической модели конвекции, включающей эндотермический фазовый переход [1-4] и процессы дифференциации мантийного вещества [5-10] (в приближении Буссинеска и изовязкой реологии). Учитывается, что: в слое D'' появляется облегченное мантийное вещество за счет отделения и перехода в ядро тяжелых металлов; в поверхностном слое экстрагируются и уходят в кору легкие составляющие, и при погружении вследствие уплотнения (эклогитизации) океанической коры образуется тяжелое мантийное вещество.
Моделирование показало [11], что эти процессы дифференциации формируют импульсно-циклический характер конвекции. Когда внизу образуется легкое вещество, а наверху – тяжелое, химические процессы вместе с термической конвекцией резко ускоряют движение, но после того, как легкое вещество всплывет, а тяжелое утонет, фактор плавучести вместе с фазовым барьером препятствует движению и способствует расслоению течения. В то же самое время латеральные химико-плотностные градиенты облегчают преодоление фазового барьера, и в численных экспериментах регулярно наблюдаются аваланши регионального масштаба. Это объясняет современную сейсмотомографическую картину мантии с прорывами 670-ти километровой границы.
Расслоенность течения ведет к охлаждению верхней мантии и нагреву нижней. В случае достижения критической (термо)плотностной стратификации в двухслойной системе один из региональных аваланшей вырастает до глобального масштаба. Другими словами, наблюдается самоорганизация интенсивной общемантийной конвекции, для описания которой хорошо подходит термин мантийный переворот [12]. Происходит организованное погружение холодного верхнемантийного вещества по одному мощному стоку, а прорывы горячего нижнемантийного вещества осуществляются в виде нескольких (3-5) восходящих суперплюмов. Такое различие объясняется эффектом заклинивания при погружении, обусловленным шаровой геометрией мантии.

Мантийные перевороты, как показывает рисунок, происходят с периодом 650-900 млн. лет (циклы Вильсона), и такой мегациклический режим конвекции определяет ступенчатый характер эволюции Земли, фиксируемый геохимиками [13]. Одностоковая структура переворотов объясняет образование суперконтинентов, их распад описывает известная модель мантийной конвекции с плавающими континентами В.П. Трубицына.
Следует отметить, что для инициации мегациклов необходимы особые условия. В наших экспериментах это было равновесное (сферически симметричное) неустойчивое по радиальной координате начальное состояние планеты. Такое начальное состояние, не только хорошо соответствует новым астрофизическим и космохимическим данным о быстрой горячей аккреции планет [14], оно также дает акцентированный стартовый мантийный переворот в форме квазикубической конвекции, что подтверждается наличием шести мощных щитов древней континентальной коры [7].
Из геологических палеореконструкций известно [7], что океаны атлантического типа формируются благодаря расколам континентов и возраст их жизни не превышает 600 млн. лет, а Пратихий океан существовал на месте современного Тихого океана с начала палеозоя. Исследование конфигураций мантийных переворотов показывает, что уже после первых мегациклов положение стока (закрывающего океаны атлантического типа и собирающего суперконтиненты) стабилизируется, а восходящие суперплюмы (создающие океаническое ложе) повторяются в противоположном полушарии. Так формируется асимметричная дипольная конфигурация планеты с пульсирующими океанами атлантического типа в одном полушарии и крупным реконструируемым на одном месте (Тихим) океаном в другом.

Литература
1. Christensen U.R., Yuen D.A. Layered convection induced by phase transitions // J. Geophys. Res. 1985. V. 90. P. 10291-10300.
2. Machetel P., Weber P. Intermittent layered convection in mantle with an endothermic phase change at 670 km // Nature. 1991. V. 350. P. 55-58.
3. Tackley P.J., Stevenson D.J., Glatzmaier G.A., Schubert G. Effects of multiple phase transitions in a 3-dimensional spherical model of convection in Earth’s mantle // J. Geophys. Res. 1994. V. 99. P. 15877–15901.
4. Трубицын В.П., Евсеев А.Н, Баранов А.А., Трубицын А.П. Мантийная конвекция с эндотермическим фазовым переходом // Физика Земли. 2007. № 12. С. 4-11.
5. Кеонджян В.П., Монин А.С. О концентрационной конвекции в земной мантии // Докл. АН СССР. 1980. Т. 253. № 1. С. 78-81.
6. Кобозев А.В., Мясников В.П. Трехмерная модель эволюции внутреннего строения планет // Докл. АН СССР. 1987. Т. 296. № 3. С. 561-565.
7. Лобковский Л.И., Никишин А.М., Хаин В.Е. Современные проблемы геотектоники и геодинамики. М.: Научный мир, 2004. 612 с.
8. Добрецов Н.Л., Кирдяшкин А.А., Кирдяшкин А.Г. Физико-химические условия на границе ядро-мантия и образование термохимических плюмов//Докл РАН 2003. Т.393. №6. С.797-801
9. Жарков В.Н. Внутреннее строение Земли и планет. М.: Наука. 1983. 416с.
10. Николаевский В.Н. Геомеханика и флюидодинамика. М.: Наука, 1996. 448 с.
11. Kotelkin V.D., Lobkovsky L.I. Numerical analysis of geodynamic evolution of the Earth based on a thermochemical model of the mantle convection: 3-D model//RJES 2004. V.6. N 6. P. 385-389 http://rjes.wdcb.ru/v06/tje04165/tje04165.htm
12. Stein M, Hofmann A. Mantle plumes and episodic crustal growth//Nature 1994. V.372. P.63-68.
13. Condie K.C. Episodic continental growth and supercontinents: a mantle avalanhce connection? // Earth Plan. Sci. Lett. 1998. V. 163. P. 97-108.
14. Kleine T., Münker C., Mezger K., Palme H. Rapid accretion and early core formation on asteroids and the terrestrial planets from Hf-W chronometry // Nature. – 2002. V. 418. – P. 952-955.

Top

Internal structure of icy satellites of Jupiter and Saturn and subsurface oceans
O.L. Kuskov Institute of Geochemistry and Analytical Chemistry RAS, Moscow, kuskov@geokhi.ru

Introduction
Data from the Galileo magnetometer produced some surprising results that suggest the presence of a conducting sub-surface layer in icy satellites of Jupiter (Europa, Ganymede and Callisto). Despite the difference in the degree of differentiation of the icy satellites, geophysical data support the possibility of the presence of a water layer with salinity comparable to seawater on Earth beneath the solid ice shell of each of them. Induced fields, observed in Europa and Callisto, indicate the strong likelihood of water oceans in these bodies.

The purpose of this study is to reproduce characteristic features of the internal structure of Europa, Ganymede and Callisto and to compare them with the internal structure of Titan. Here, taking into account ordinary and carbonaceous chondrites as representatives of nebula matter and the reliable estimates of moment-of-inertia factors (MOI), we also concentrate of assessing various geochemical parameters of the satellites. The phase diagram of H2O and equations of state of high-pressure ices and meteoritic matter are taken into account. The general methodology is to combine the geophysical and geochemical constraints and thermodynamic approach, and to develop, on this joint basis, the self-consistent models of Europa, Ganymede and Callisto, accounting for their composition and internal structure [1].

The model and geophysical constraints
The problem of modeling the internal structure of a satellite is described by a system of equations specifying the conditions of thermodynamic and hydrostatic equilibrium, equations of state of phases, equations of heat transfer by heat conduction, and mass and moment conservation. The mass, average density and moment of inertia are used as input data for determination of (1) the thickness and phase state of an outer water-ice shell, (2) the density distribution with depth, and (3) the core sizes and masses. The density variations in the mantle and core radii are found by the Monte-Carlo method.

Geochemical constraints
Additional geochemical limits on mantle density (and indirectly on the temperature) are introduced from the density of the equilibrium phase assemblages (phase diagrams) calculated from the composition of silicate fraction of ordinary (H, L and LL) and carbonaceous (minus volatiles) chondrites. The mantle densities varied within the density continuum in the range of 3320-3670 kg/ m3, encompassing the possible density variations of mineral assemblages of the silicate fraction of chondrites at 0.2-4 GPa and 300-1400ºC.

Thermodynamic constraints
The phase compositions and mantle densities are modeled within the system Na2O-TiO2-CaO-FeO-MgO-Al2O3-SiO2 including the solid solution phases (olivine, spinel, ilmenite, garnet, orthopyroxene, and clinopyroxene). The equilibrium phase assemblages were calculated using the technique of free energy minimization. The equations of state of minerals were calculated in the Mie-Grüneisen-Debye approximation.

Model of a core Three boundary models are considered for a satellite core: a Fe-10 wt%S, a eutectic Fe-FeS and a FeS core. The density of 5700 kg/m3 is adopted for an Fe-10 wt%S alloy at 5 GPa and 1500oC; ρ=4700 and ρ5150 kg/m3 are adopted for a central FeS core and for a eutectic Fe-FeS.

Water-ice shell
The densities of liquid and high pressure ices were calculated by the equations of state in accordance with the phase diagram of H2O.

Europa
The results show that Europa is differentiated into a water-ice shell, anhydrous mantle and iron-sulfide core. Both L/LL and CM chondrite compositions match the total mass and moment of inertia value of Europa and can be regarded either as the primary material of Europa (carbonaceous chondrites) or as a reasonable analogue for its anhydrous rock-iron core (ordinary chondrites). However the amounts of metallic iron and iron sulfide, and the Fetot/Si ratio of Europa's anhydrous rock-iron core are not consistent with the bulk compositions of the most oxidized CI chondrites and the most reduced H chondrites. It is likely that Europa inherited a significantly higher proportion of material close to the moderately oxidized L/LL type chondrites rather than to the carbonaceous chondrites. The allowed thickness of Europa's H2O layer (whether liquid or ice) ranges from 11510 km (6.80.6% of total mass) for a differentiated L/LL-type chondritic mantle with a crust to 13510 km (7.90.5%) for an undifferentiated L/LL chondritic mantle.

Ganymede
Two alternative density models of an outer shell are considered. Model (A) - an outer shell is completely composed of the high-pressure ice phases (no water is present), resulting in a maximum in the density of an outer shell. Model (B) - in the three-layer model of an outer water-ice shell, we assume that below a shell of ice-I (30-120 km thick), a liquid layer of 230-140 km thick may exist, resulting in a minimum density of an outer shell. The ice-V + ice-VI + water triple point lies at 273 K and 0.64 GPa. We adopted a “conductive” model where a mixed layer of water and high-pressure polymorphs of ice may coexist at depths between 260 km and an ice-rock interface. Our calculations show that the ice thickness of the outer shell in model (A) is about 890-920 km and in model (B) is 780-850 km.

Callisto
We show that Callisto must only be partially differentiated into an outer ice-I layer, a water ocean, a rock-ice mantle, and a rock-iron core free of ice (mixture of anhydrous silicates and/or hydrous silicates + Fe-FeS alloy, 3150 <  < 3620 kg/m3). Assuming conductive heat transfer through the ice-I crust, heat flows were estimated and the possibility of the existence of a water ocean in Callisto was evaluated. The liquid phase is stable (not freezing) beneath the ice crust, if the heat flow is between 3.3 and 3.7 mW/m2, which corresponds to the heat flow form radiogenic sources. The thickness of the ice-I crust is 135-150 km, and that of the underlying water layer, 120-180 km. The allowed total (maximum) thickness of the outer water-ice shell is up to 270-315 km. The results of modeling support the hypothesis that Callisto may have an internal liquid-water ocean.

Conclusion
(1) Models of the internal structure of completely differentiated Europa and Ganymede, and partially differentiated Callisto have been constructed on the basis of Galileo gravity measurements, geochemical constraints on composition of ordinary and carbonaceous chondrites, and thermodynamic data on the equations of state of water, high-pressure ices, and meteoritic material. Comparison of the internal structure of the Galilean satellites and Titan has been made [1].

(2) Geophysically and geochemically permissible thicknesses of Europa's and Ganymede’s outer water-ice shell are determined. We show that Callisto must only be partially differentiated into an outer ice-I layer, a water ocean, a rock-ice mantle, and a rock-iron core. The results of modelling support the hypothesis that Callisto may have an internal liquid-water ocean.

(3) The correspondence between the density and moment of inertia values for bulk ice-free Io, rock-iron core of ice-poor Europa, and rock-iron cores of Ganymede and Callisto shows that their bulk compositions may be, in general, similar and may be described by the isochemical composition close to a material of the L/LL type chondrites. Therefore, planetesimals composed of these types of ordinary chondrites could be considered as analogues of building material for the rock-iron cores of the Galilean satellites. Similarity of bulk composition of the rock-iron cores of the inner and outer satellites implies the absence of iron-silicon fractionation in the protojovian nebula.

1. Kuskov O.L., Dorofeeva V.A., Kronrod V.A., Makalkin A.B. (2009) The systems of Jupiter and Saturn: Formation, composition and internal structure of large satellites. Мoscow, LKI, 576 pp.

Top

Constraining the composition and thermal state of the cratonic mantle from inversion of seismic data
O.L. Kuskov Vernadsky Institute of Geochemistry and Analytical Chemistry RAS, Moscow, kuskov@geokhi.ru

Introduction
Quantitative estimation of the temperature distribution in the Earth’s mantle is a key problem in petrology and geophysics. Geophysical and geochemical data alone are insufficient to constrain mantle temperatures. Seismic and thermal studies show both seismic heterogeneity and substantial variations in the distribution of temperature in the continental upper mantle. Thermal models are uncertain because of the large uncertainties in the values of crustal heat production. Despite the uncertainties in geotherms predicted from seismic models, seismic velocities may provide an alternative source of information of thermal regime of the upper mantle.

In this study, we discuss the method of estimating temperature, composition and thickness of the lithosphere beneath the Kaapvaal and Siberian cratons from seismic velocities. The thermal boundary of the lithosphere is defined by the intersection of the calculated temperature profile in the conductive region with the potential adiabat. Based on self-consistent thermodynamic approach, we infer the temperature distribution models at 100 to 300 km depth for the upper mantle beneath the Kaapvaal and Siberian cratons from absolute P- and S-wave velocities and geochemical constraints (garnet-bearing lherzolite xenoliths, average composition of garnet peridotite and primitive mantle composition.

Input data and method of solution
The phase composition and physical properties of the lithospheric mantle were modelled within the Na2O-TiO2-CaO-FeO-MgO-Al2O3-SiO2 (NaTiCFMAS) system including the non-ideal solid solution phases: olivine (Ol), spinel (Sp), plagioclase (Pl) and ilmenite (Ilm) – binary solutions; garnet (Gar: almandine, pyrope, grossular); orthopyroxene (Opx: MgSiO3, FeSiO3, Ca0.5Mg0.5SiO3, Ca0.5Fe0.5SiO3, Al2O3); clinopyroxene (Cpx: same components as in Opx plus jadeite end-member). Chemical reactions in this system are independent of oxygen fugacity. We include TiO2 in the ilmenite phase but exclude some minor elements as less significant for seismic properties of the deep interior. By applying suitable thermodynamic models for the equation of state of minerals and solid solutions and the technique of free energy minimization, phase relations, seismic velocities and density can be derived entirely from an internally consistent model, with only thermodynamic data and bulk composition as input. The bulk composition (unless noted otherwise) is assumed fixed with depth whereas the phase composition, phase proportions and chemical composition of the coexisting phases are changing with pressure and temperature. Input data for the self-consistent thermodynamic quantities (enthalpy of formation, entropy, heat capacity, and Margules parameters for non-ideal solid solutions), for the EOS of minerals as well as for shear moduli and their derivatives are summarized in our THERMOSEISM data base, which includes a series of computer codes for both forward and inverse solutions of the problem of mantle structure. Temperature dependence of seismic velocities comes both from anharmonicity and anelasticity; the latter is an important effect at seismic frequencies. Following many authors, anelastic effects were included for comparing theoretical velocities with geophysical interpretations. Forward code calculates the equilibrium phase assemblages and their velocities and densities. Inverse code computes the temperature distribution in the mantle from seismic and compositional constraints. Anharmonic and anelastic effects are taken into account. The output results contain the self-consistent information on phase assemblages, densities and velocities. The approach used here requires a small number of thermodynamically defined parameters and has important advantages over earlier procedures, which contain no information about entropy, enthalpy and Grüneisen parameter [1].

Kaapvaal and Siberian cratons
The great majority of mantle xenoliths are peridotites. Olivine and orthopyroxene are the main mineral constituents of these rocks. Other lithologies that occur in the upper mantle are significantly less abundant than peridotites. Here we restrict our analysis to some of these xenoliths and use also the average composition of garnet peridotite (GP model). We approximate the composition of the “asthenosphere” by primitive mantle composition (PM model), which is close to the pyrolite composition. GP compositions have higher Mg numbers and are low in Al and Ca, and garnet and clinopyroxene compared to the primitive mantle.
We have used two compositional models for the lithospheric mantle beneath the Kaapvaal and Siberian cratons: a fixed composition of garnet peridotite (depleted in basaltic components) and variable composition from the GP model to the PM fertile model at depths of 180-250 km. We have used the fertile composition of primitive mantle for the normal mantle

Calculation of density and velocities from petrologic models (forward modelling)
P- and S-velocities and densities as a function of depth are calculated for the depleted garnet peridotite (GP model) and non-depleted primitive mantle composition (PM model), and for the garnet lherzolite xenolith samples from their bulk chemical compositions at equilibration pressure and temperature.

Inferring upper-mantle temperatures from seismic and geochemical constraints (inverse modelling)
We inverted for temperature the recent P and S velocity models of the Kaapvaal craton as well as the IASP91 reference Earth model. Several long-range seismic profiles were carried out in Russia with Peaceful Nuclear Explosions (PNE). The velocity models from PNEs recorded along these profiles were used to infer upper mantle temperature profiles beneath the Siberian craton. The seismic profiles were inverted on the basis of low and high temperature xenoliths of garnet peridotites from kimberlite pipes of the cratons in order to gain insights into the temperature sensitivity to variations in the composition and mineralogy of xenoliths. Such a test can provide constraints on the compositional (vertical and lateral) heterogeneity of the upper mantle.

Results
Our results show that chemically depleted lithosphere beneath the Kaapvaal craton exists to depths of about 175-200 km. At greater depths (~200-250 km), the lithospheric material appears to be substantially enriched in FeO and refractories (Al2O3 and CaO) compared to the depleted GP compositions (i.e., low-T and high-T xenoliths) but simultaneously also must be more depleted than primitive mantle. Thus, the mantle composition in a depth range of 200 and 250 km is more enriched than the overlying mantle but it is more depleted than underlying fertile mantle. We arrive to a conclusion about the impossibility of a uniform composition throughout the entire Kaapvaal cratonic mantle. This may suggest a strong compositional gradient below ~200 km and a vertical zonation in the composition of the Kaapvaal lithospheric mantle. A lithospheric mantle would be expected to be chemically stratified, becoming more fertile in composition with depth. At depths of about 275 km, Kaapvaal cratonic mantle does not differ from normal mantle. We thus derive a lithosphere thickness of roughly 275 km for the Kaapvaal craton [1, 2].

2D thermal profiles of the Siberian craton inferred from the available seismic profiles are discussed. The boundary of the thermal lithosphere approximately conforms to 1450oC isotherm that corresponds with depth of 300-350 km.

References
1. Kuskov O.L., Kronrod V.A., Annersten H. (2006) Inferring upper-mantle temperatures from seismic and geochemical constraints: Implications for Kaapvaal craton. Earth Planet. Sci. Lett., 244, 133-154.
2. Kronrod V. A., Kuskov O. L. (2007) Modeling of the thermal structure of continental lithosphere. Izv. Phys. Solid Earth, 43, 91–101.

Top

The effect of climate and damage on generating plate tectonics on a planet
William Landuyt^1,2, David Bercovici²
1. IGPP, Scripps Institution of Oceanography, UCSD, San Diego, CA USA
2. Department of Geology and Geophysics, Yale University, New Haven, CT USA

The generation of plate tectonics on Earth and its absence on the other terrestrial planets (especially Venus) remains a significant conundrum in geophysics. We propose a model for the generation of plate tectonics that suggests an important interaction between a planet’s climate and its lithospheric damage behavior; and thus provides a simple explanation for the tectonic difference between Earth and Venus. We propose that high surface temperatures will lead to higher healing rates (e.g. grain growth) in the lithosphere that will act to suppress localization, plate boundary formation, and subduction. This leads to episodic or stagnant lid convection on Venus because of its hotter climate. In contrast, Earth’s cooler climate promotes damage and plate boundary formation. The damage rheology presented here attempts to describe the evolution of grain size by allowing for grain reduction via deformational work input and grain growth via surface tension-driven coarsening. We present results from convection simulations and a simple “drip-instability” model to test our hypothesis. The results suggest the feasibility of our proposed hypothesis that the influence of climate on a lithosphere with a damage rheology may control the mode of tectonics on a planet.

Top

Comparison between mantle seismic tomography and past plate boundaries of the western Pacific
So, Byung-Dal (qudekf1@snu.ac.kr) and Lee, Sang-Mook (smlee@snu.ac.kr) Seoul National University, School of Earth and Environmental Sciences, Seoul 151-747, Korea

The western Pacific is and has been host to many important plate boundaries. These include subduction zones and back-arc basins many of which have not only appeared and disappeared but also have changed their positions and dimensions. We know from the present day bathymetry and seismicity where the active boundaries are located, but it is not straightforward to determine where the boundaries have been in the past. In the case of subduction zones, the difficulty is compounded by the fact that, unlike mid-ocean ridges or marginal basins, they do not leave geological traces on the seafloor. Yet understanding both the present and past plate boundaries are important for understanding evolution of western Pacific. In this study, we use a combination of existing plate kinematic models and seismic tomography data to identify past subduction zones and their evolution. For instance, it is unclear whether the divergent plated boundary between Caroline plate and Philippine Sea plate also known as Ayu Trough was the subduction zone consuming the westward migration of Pacific plate in the past or if this movement was accommodated by subduction zones along the eastern boundaries of Caroline plate. Although past subduction zones might have disappeared from the seafloor, they may still leave important traces such as velocity anomaly in the mantle, which can be tracked back in time. By carefully comparing the presumed location of past plate boundaries with seismic tomography images, we intend to examine the validity of different plate kinematic models and improve our understanding of the geological evolution of the western Pacific.

Top

Thermodynamics of Earth Surface Instabilities
Klaus Regenauer-Lieb(1),(2), David A. Yuen(3), Florian Fusseis(1) email: klaus.regenauerlieb@csiro.au, daveyuen@gmail.com, fusseis@cyllene.uwa.edu.au
1. Multi-Scale Earth System Dynamics, School of Earth and Environment, University of Western Australia, Crawley, WA, 6009 Australia
2. CSIRO Exploration and Mining, PO Box 1130, Bentley WA 6102, Australia, FAX ++61864368555*
3. Dept. of Geology and Geophysics and Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, MN,55455 U.S.A

The total rate of rock deformation results from competing deformation processes, including ductile and brittle mechanisms. Particular deformation styles arise from the dominance of certain mechanisms over others at different ambient conditions. Surprisingly, rates of deformation in naturally deformed rocks are found to cluster around two extremes, representing co-seismic slip rates or viscous creep rates. Classical rock mechanics is traditionally used to interpret these instabilities. These approaches consider the principle of conservation of energy. We propose to go one step further and introduce a nonlinear far-from-equilibrium thermodynamic approach in which the central and explicit role of entropy controls instabilities.

Classical Mechanics
Our approach is based on a classical mechanics framework, which is cast here in a consistent thermodynamic description. The new extensions to this theory allow us to explore time-dependent thermodynamic feedback processes which rely on consideration of the explicit role of entropy production described in section entitled “Time-Dependent Thermodynamics”.

Rheology
We describe in the following the classical reasoning for using the additive strain-rate decomposition of elasto-visco-plastic rheology. In geodynamics the motivation for this serial implementation does not stem from a thermodynamic perspective. However, it is the most practical approach from a thermo-mechanical perspective because it separates the Helmholtz free energy from the dissipation potential. The Helmholtz free energy is a path-independent functional of the state variables such as temperature, strain rate and pressure, while the dissipation potential corresponds to the path-dependent entropy production (Collins, 2005).

The first term in the bracket on the right hand side describes the elastic strain, where E is Young’s modulus, ν is Poisson ratio and λ is the coefficient of thermal expansion.

is the objective co-rotational stress rate and

is the Kronecker delta. The second term on the right hand side describes the plastic strain, which in our simple formulation scales linearly with the trace of the Cauchy stress tensor

. The last term on the right hand side describes the creep strain where A and n are power-law material constants, Q is the activation enthalpy and R is the universal gas constant. J₂ is defined as the second invariant of the deviatoric stress tensor

and

is the deviatoric stress

.

The rheology used here is based on a classical mechanical approach. Thermodynamic consistency is ensured by consideration of Ziegler’s maximum dissipation principle (Ziegler, 1983), which postulates that there is no work-free plastic deformation. Mathematically this translates into the postulate of orthogonality between the dissipative stress and the dissipation potential.

Conservation of mass
The conservation of mass is also known as the continuity equation:

where

is the density and u is the local material velocity vector of the volume under consideration; the second term describes the divergence of the velocity field. The density-time derivative is implicitly encapsulated in an isentropic thermodynamic formulation of the energy balance law. This is done here in terms of an equation of state, however, in geodynamics the density is often considered to be independent of time.

Conservation of momentum
We also require the classical momentum equation

where

is the divergence of the Cauchy stress tensor and f is the body force.

Energy Conservation
The classical linear irreversible thermodynamic formulation in the case of thermomechanical, creeping flow or quasi-static equilibrium is written in terms of the first law of thermodynamics as the energy balance

Where u is the specific internal energy, q is a heat flux vector, n is a normal vector to the surface

of the volume (V) and

is the rate of total work. We use the “tilde” for the work rate to emphasize that time integrals are path dependent. The term r_i collects all internal heat sources such as radiogenic heat production.

Entropy Conservation
The internal irreversible entropy production is given by:

where

is the path dependent total specific entropy change rate in the local form,

_irr is the irreversible specific entropy production rate, and T the absolute temperature. Through the introduction of the irreversible entropy production equation (5) we make use of the concept of writing the second law of thermodynamics in an equality form (Tolman and Fine, 1948). It is obvious from equation (5) that it is the entropy rather than the temperature that describes the underlying physical processes. The first term on the right hand side is the heat source through the entropy production caused by mechanical work. This equation shows the important link of mechanical deformation and heat production, which needs to be considered in a fully coupled way. Classical mechanics decouples the heat problem from the mechanical problem, and therefore has to use empirical observations (weakening laws, localization mechanisms etc.) from experiments to derive time dependent processes. In geodynamics the time frame is too long for making reproducible experiments and we have to consider the explicit time-dependent coupling of heat supply, rate of mechanical work and rate of internal energy supply.

Time-Dependent Thermodynamics
In addition to the above-described classical mechanical methodology, we explore the explicit role of entropy production. The entropy production is related to the rate of change of internal energy through the Helmholtz formulation:

where

is the specific Helmholtz free energy dependent on the state variables. We differentiate equation (6) with respect to time and obtain:

The second term on the right hand side describes the rate of entropy production and the third term describes the coupled variation of temperature with time. Feedback between the three terms on the right hand side of equation (7) causes localization phenomena in the models, such as shear bands, stemming from feedback effects and leading to self-regulation.

We define the strength measure of the lithosphere as the time integral of the total recoverable power

that can be stored in the thermodynamic system, where a_ijdescribes all recoverable state variables, like elastic strain. The weak (local)form of this power is given by:

The strong form of this equation can be found in (Regenauer-Lieb et al., 2009). The thermal expression of the irreversible entropy production is described by

where

is the feedback term due to thermal-mechanical coupling. Here we consider thermal-elasticity where

is the thermal-elastic temperature change. This term incorporates the thermal expansion coefficient λij which we choose to be isotropic . Other possible applications of this term are the latent heat release assuming that the state variable is the fractional volume of a given phase. The irreversible entropy production is (Regenauer-Lieb et al., 2009):

The first term on the right is the total rate of mechanical work, the second term is the fraction of the work that is stored - in our implementation the elastic work rate - and not converted into heat. In a more general formulation for the state variables this term could be the power stored through chemical changes, defects and other microstructures. Combining equations (9) and (10) and rearranging for the rate of temperature change, we obtain the energy equation with additional thermalmechanical coupling terms (first three terms on the right hand side of the equation):

These explicit time dependent terms are the sources of the instabilities causing the strain rate cascade of surface instabilities from slow geodynamic rates to slow earthquakes.

Time-dependent Earth Surface instabilities
Observational evidence for the different mechanisms come from a wide range of disciplines, from extrapolation of experimental data, micro-structure, microgeochemistry (isotope measurements, 'geospeedo-metry'), remote sensing(GPS/InSAR) and geophysical measurements (borehole strainmeters). Total rates of natural deformation processes have been bracketed (Figure 1) either at seismic rates (e.g. fracture propagation, frictional sliding/melting, which happens within seconds) or at strain rates of 10-10 to 10-15 s-1 where point/line-defect diffusion controls the viscous (i.e. temperature-activated) response of a rock. A geological system where deformation takes place at all strain rates is the seismogenic crust. Deformation occurs on the time-scale of the earthquake cycle on the order from tens to hundreds of years. During interseismic periods various ‘slow’ diffusion processes ultimately control most deformation, whereas during earthquakes ‘rapid’, unstable brittle processes control the deformational response of the fault to co-seismic stress release. Consequently, during the earthquake cycle, deformation rates vary over about 14 orders of magnitude. Within the earthquake cycle "slow earthquakes" release strain energy over a time scale of several hours to days, often as a precursor to larger earthquakes; e.g. [Linde et al., 1996]. This time-scale suggests that slow earthquakes happen at intermediate rates between those achieved during inter-seismic and coseismic periods.

Figure 1 Strain rates symbols: GPS (gray), InSAR(white), Strainmeters (black). Triangles show postseismic deformation from the Landers earthquake, squares show slow earthquakes, open squares single station observations. GPS and InSAR data have higher sensitivity at long periods (months to decades) while borehole strain meters are better at short periods. While some data is derived close to their detection limit, the plots illustrates the point that there is no separation between the classical geodynamic strain rates (dashed line plate motion) and seismological strain rate. Geological deformation happens at all rates. SOURCE: PBO Steering Committee, The Plate Boundary Observatory: Creating a four-dimensional image of the deformation of western North America, White paper providing the scientific rationale and deployment strategy for a Plate Boundary Observatory based on a workshop held October 3-5, 1999

This data shows that there is a clear need to fill the gap between geodynamic and seismological modeling. In Figure 2 we compare the time dependent processes of the thermodynamic approach arising out of shear-heating and thermal-expansion feedback with the classical quasistatic approach which is independent of time. The geodynamic simulation based on the thermodynamic approach reproduces unstable slip range from 1000 years to million of years as suggested by Figure 2.

Figure 2 Two examples for a crust in extension: a) Time dependent thermodynamic approach and its creep instabilities b) compared to the classical mechanics approach courtesy of Chris Wijns. The thermodynamic approach shows a much richer patterns in strain a) and strain rate b). The magnitude of the slip increases in a power law fashion for shorter time periods.

References
Collins, I.F., 2005. Elastic/plastic models for soils and sands. International Journal of Mechanical Sciences, 47(4-5 SPEC. ISS.): 493-508.

Regenauer-Lieb, K., Yuen, D. and Fusseis, F., 2009. Landslides, ice quakes, earthquakes : a thermodynamic approach to surface instabilities. Pure and Applied Geophysics:166, 1-24.

Tolman, R. and Fine, P., 1948. On the irreversible production of entropy. Reviews of Modern Physics, 20: 51-77.

Ziegler, H., 1983. An Introduction to Thermomechanics. North Holland, Amsterdam.

Top

Deformation and seismicity associated with rift zones crossing passive continental margins
Lyakhovsky V. ⁽¹⁾, Segev A. ⁽¹⁾, Weinberger R. ⁽¹⁾, Schattner U. ⁽²⁾
(1) Geological Survey of Israel, 30 Malkhe Israel St., Jerusalem 95501
(2) Straus Department of Marine Geosciences, CSMS, Faculty of Science and Science Education, University of Haifa, Haifa 31905

We present a 3-D numerical modeling at lithosphere scale to assess the impact of the transition from the thick continental crust to the thin oceanic crust along passive continental margins on the style of rift propagation. During the last decades much attention has been paid to the modeling of the extensional tectonics and rifting of continental and oceanic lithosphere. Vink et al [1984] noted that thick continental lithosphere, with its granitic crust, should be weaker than oceanic lithosphere. They also noted the role of the depth to Moho and strength of the upper mantle, which is dominated by the temperature-dependent olivine rheology. Thus, a rift might stop propagating toward oceanic lithosphere at a continental margin and might switch tectonic activity. This has been suggested for the NW continuation of the active Red Sea-Suez rift system and initiation of the Dead Sea Transform (Steckler and ten Brink, 1986). With the onset of the Red Sea opening (about Oligocene times) the sub-parallel Azraq-Sirhan rift was also activated and propagated in a NW direction from the Arabian plate toward the Levant oceanic basin. The gradual thinning of the crystalline crust and strengthening of the lithosphere in the direction of the Azraq-Sirhan system along the Lower Galilee might lead to a decrease in the rate of rift propagation, or even terminate the rifting above the Levant continental margin.

We address these questions utilizing 3-D lithosphere-scale numerical simulations of the plate motion and faulting process base on the damage mechanics. The 3-D modeled volume consists of three main lithospheric layers: an upper layer of weak sediments, middle layer of crystalline crust and lower layer of the lithosphere mantle. The total strain in each layer is the sum of (1) elastic strain, (2) damage-related inelastic strain, and (3) ductile strain. We use a visco-elastic damage rheology model to calculate the elastic strain coupled with evolving material damage and damage-related inelastic strain accumulation. Newtonian viscosity governs the ductile strain component in the sedimentary layer, while temperature-dependent power-law rheology is used for the ductile strain component in the lower crust and upper mantle. The variables controlling the kinetics of the damage evolution are constrained by results of laboratory analyses of stress-strain relations, acoustic emission and frictional data, as well as by scaling relations between seismic moments and rupture areas. The material constants for a creep law of a diabase are used to control the ductile flow of the lower crust and a ductile flow of a wet olivine for the upper mantle below the Moho interface. Depth-dependent temperature distribution in the simulated volume follows a given geothermal gradient corresponding to the regional heat flux. Stresses in the modeled volume are dictated by lithostatic pressure, mass density, elastic moduli and tectonic loading due to a plate motion far from the simulated region. Neither constant velocity nor constant force boundary conditions are appropriate for long-term simulations which aim at studying the evolution of earthquakes and faults over many large earthquake cycles. The boundary conditions are updated during the simulations to account for the evolution of the elastic properties within a simulated volume, as well as finite tectonic stresses which correspond to the plate motion. To perform long-term simulations with appropriate boundary conditions at the edges of a finite-size region, we use variable boundary forces proportional to the stiffness of virtual springs multiplied by the mismatch (slip-deficit) between the far field plate motion and the displacement of boundary nodes.

Results of the modeling demonstrate how the depth to the Moho interface affects the geometry of the propagating rift system and its associated strain field as well as its seismicity pattern. In the case of a homogeneous lithosphere with a flat layered structure, the rift system bisects the modeled area in a direction corresponding to the remote loading. With the same boundary conditions and physical properties of rocks as in the case of the flat structure, a rift terminates above the passive continental margin and a new fault system is created normal to the direction of the rift propagation. These results demonstrate that the local lithosphere structure is one of the major factors controlling the geometry of the evolving rift system, faulting and seismicity pattern.

For the presented analyses of the Azraq-Sirhan rift system the lithosphere structure of the study area was reconstructed to the end of the Eocene time. The numerical results demonstrate a NW propagation of the rift due to remote loading applied perpendicular to the rift direction. A new NE-trending fault system developed and the rift terminated next to the continental passive margin. Numeric simulations were performed with various prescribed regional heat fluxes, ranging between 40 mW m-2 and 60 mW m-2, corresponding to cold and normal lithospheres. In the cold lithosphere the localization of the propagating rift zone was fast and terminated above a thicker crust (ca. 30 km), close to the present location of the Dead Sea transform. In the hotter lithosphere, the rift propagated further northwestwards toward the oceanic lithosphere and terminated across a thinner crust. The pattern of the fault systems close to the northwestern tip of the propagating rift, demonstrates a fan shaped structure.

Top

ON THE POSSIBILITY OF CRATER FORMATION ASSOCIATED WITH AN ASCENDING PLUME
A.B. Medvedev Russian Federal Nuclear Center– VNIIEF, Sarov, root@gdd.vniief.ru

One of the peculiar features (fig.) of experimental data about shock compression of silicates consists in the fact that in the coordinates of pressure Р – specific volume V shock adiabats of a substance having a decreased initial density (the initially porous samples characterized by the initial temperature 300 K, or, in case of fayalite, the initially liquid samples at 1573 K) are located far more to the left of the shock adiabat of a substance conforming to higher values of (solid samples with an initial density close to the crystal density and 300 K) in a wide pressure range. Ordinary equations of state (EOSs) of condensed material can not describe a similar peculiarity within a wide range of states, even if EOSs take into account phase transitions in the crystalline state and melting (because of their relative locality). This unusual property may be qualitatively explained through amorphization of silicates at the increased Р,Т-conditions and through the fact that an amorphous substance (glasslike and liquid) has a wide range of states characterized by a negative value of a coefficient of thermal expansion (at Р=const. a temperature of a shock-compressed initially porous substance is higher than a temperature of a shock-compressed solid substance, therefore an abnormal relative position of adiabats suggests that ). In this case, there is an analogy with an anomaly for water in the range of t 0-4 оС at a qualitative level. The anomaly gives way to a normal ( ) behavior with increase in Т,Р- conditions. For this reason, in particular, the shock adiabats SiO2 =1.35, 0.8, 0.55 g/cm3 at maximum Т,Р- parameters are placed (fig.) to the right of the adiabat of fused quartz =2.204 g/cm3. Obviously, the possible existence of such unusual properties of silicates should be taken into account in modeling processes in the mantle of planets (what is very significant for material buoyancy), at least, hypothetically. Such an approach was applied in [5] at a qualitative level with the use of a model EOS of (the shock adiabats calculated on the basis of EOS are presented through lines in fig.) for Mercury and Moon in the framework of the hypothesis about the existence in the distant past in their interior (at the mantle-core boundary) of sources of powerful thermal plumes. The use of this EOS to their mantle points to the possibility of anomality of average states of a mantle material in the range of Т 1500-3000 K, Р 0-10 GPa and normalization of states of a overheated substance of thermal plumes. Under similar conditions the anomaly (a) restricts the buoyancy of plumes (a hot material must be overheated sufficiently for ascending, for the purpose of overcoming a density barrier caused by anomaly, in this case, the ability to the ascent increases with increasing a column height of a hot substance; as a result, only overheated columns comparable in height with the vertical dimension of the mantle, can ascend from a bottom to a crust along channels at moderate overheating 1000 K), (b) enables the ejection (pushing out) of a cold peripheral substance (surrounding an emerged hot column) from a medium, (c) causes compaction (under horizontal heat exchange) of a vertical contact zone of a column and surroundings at the first temporary stages after its ascent as well as (d) results in depth compaction of the mantle due to the long-term heat supply. Such phenomena can lead to vertical craterlike deformations of the crust in areas of ascending large plumes. The above-mentioned feature “a” is analogous to the restricted buoyancy of hot water in cold water, characterized by t=0 оC (only the hot element of water characterized by t>8 оС can emerge in cold water), and properties “c”, “d” are due to the fact that when cooling the emerged hot element up to a warm state (characterized by t 0-8 оC) it becomes denser than surrounding cold water, characterized by t=0 оC. As a result of this, a level of warm water is reduced (in a tube, whose wall play a role of a limiter for spreading of the emerged element) with relation to a level of cold water, although immediately after the emergence of the hot t>8 оC element this level was higher or equal to it. Significant implications of such an anomaly for geophysical processes can also be possible. In particular, if the substance of the upper mantle is characterized by the average , and hot columns-plumes characterized by , comparable in size to the upper mantle, penetrate into it from the lower mantle, the substance will be compacted (at horizontal heat exchange) at the contact zone between the hot column and the surrounding upper mantle after their emergence to the crust (both hot and cold materials and their mixture are compacted after heat is transferred from the hot material into the cold material), and the column will be compacter than a medium after cooling up to a warm state (as in the case of water). If mountains (or ocean islands) and rifts are considered as a consequence of thermal columns-plumes emerged along a linear chain, the mentioned compaction (controlled by means of the velocity of the horizontal heat exchange) of the area (areas) of the contact zone can bring about a depth seismicity (step-wise movements of non-ideally flowable material inside and in the vicinity of the area of compaction) of similar regions. The velocity of the wave of the horizontal heat exchange over a interval of t~1 Myr is estimated at ~1 cm/yr from the relation (where x– the distance traveled by the heat wave, t- the time after the plume emergence, ~0.01 cm2s – the thermal diffusivity). References: (1) Chen G.Q. et. al. Physics of Earth and Planetary Interiors. 2003. V.134. р.35; (2) LASL Shock Hugoniot Data. Ed.S.P.Marsh. Berkeley-Los Angeles-London: Univ. California Press.1980; (3) Jackson I., Ahrens T.J. J. Geophys. Res. 1979. V.84. №B6. р.3039; (4) Simakov G.V., Trunin R.F. Izvestiya. Akad. Nauk SSSR. Fiz. Zemli. 1990. №11. p.72; (5) Medvedev A.B. Fizika Zemli. 2008. №4. p.48.

Top

Generation of intermediate to deep earthquakes by self-localizing thermal runaway: insights from petrological and numerical studies
Sergei Medvedev¹, Timm John1², Lars H. Rüpke1³, Yuri Podladchikov¹, Torgeir B. Andersen¹, Hakon Austrheim¹
1. Physics of Geological Processes PGP, University of Oslo, Norway
2. Institut für Mineralogie, University of Münster, Germany
3. The Future Ocean, IFM-GEOMAR, University of Kiel, Germany

Convergent margins are characterized by a strong seismic activity with earthquakes occurring at depths of up to 700km. Shallow earthquakes (<50km) are explainable by the brittle failure of rocks. At greater depths, however, the increased ambient pressure should inhibit brittle failure. We present the results of a joint petrological, analytical and numerical study of thermal runaway as a mechanism for intermediate to deep earthquakes (John et al., 2009). The Precambrian Kråkenes gabbro, located within the Caledonian (410±10Ma) high- to ultrahigh pressure metamorphic Western Gneiss Complex (WGC) in Norway, provides a spectacular example of intensely localized shear failure of rocks at high confining pressures. The key field observations are: (i) co-facial, eclogite-facies shear zones and pseudotachylytes coexist and share the same structural orientation; (ii) both have higher degrees of hydration, caused by infiltration of external fluids, and up to 3-orders of magnitude smaller grain sizes than the almost dry wall rock; and (iii) both types formed in a single, continuous and fast event. Using analytical and numerical models we explore under which conditions (e.g. differential stress, P-T, and rheology) these shear zones/pseudotachylytes may have formed and why they coexist within the same structural orientation. We demonstrate that the possible failure mechanism is thermal runaway of ductile shear zones within visco-elastic slabs. Irreversible ductile deformation dissipates mechanical work into heat that can lead to thermal softening and failure by progressively self-localizing thermal runaway (SLTR, Braeck and Podladchikov, 2007, Braeck et al., submitted)

Two numerical simulations were designed to analyze how gabbro reacts to imposed stress and to demonstrate that a one-dimensional model of thermal-runaway is consistent with the field observations. Both simulations used an ambient temperature of 680ºC to be in agreement with the ambient metamorphic PT conditions, an initial differential stress of 1.5 GPa, and the rheology of diabase (which is rheologicaly equivalent to a gabbro in applied conditions). The models share the same geometry with a perturbed zone much smaller than the model domain, in accordance with the field observations. The only difference between the two models is in the slightly lower (<1%) initial effective viscosity of the perturbed zone in simulation (2). Figure 1 shows the simulated deformation pattern at four different times. During the first phase (up to time C) the two simulations develop similarly and create similar shear zones. The transition from time C to D in simulation (2), results in catastrophic thermal-runaway with localization (>103 in the central part of the zone, which is <1 mm wide), a stress drop of ~750 MPa in milliseconds, extensive shear heating that melts the deformation zone (D2), and no subsequent deformation. The thermal field in simulation 2 continues to evolve after the initial melting (D2-E2), and this is consistent with the observed microstructures that indicate stress free crystallization. The numerical simulations provide strong support for the idea that SLTR is a viable mechanism for intermediate-depth earthquakes, under conditions corresponding to those of the metamorphic assemblages in the studied gabbro. An initial stress of 1.5 GPa is chosen in our simplest models. More elaborate simulations that include continuous stress loading and weakening due to fluid infiltration or low temperature plasticity at high confining pressures could significantly lower the stress required to initiate SLTR.

Our simulations provide an explanation for the observed coexistence of ductile deformation (shear zones) and brittle-like catastrophic failure (pseudotachylyte veins) with seismogenic strain rates and active melting at the deep-earth conditions. We demonstrate that small rheological differences in the system determine whether a shear zone will develop at a non-seismic velocity or if it will eventually rupture seismically. Our findings imply that at greater depths (>60-80 km), failure by SLTR is more likely than brittle failure, and that SLTR is a viable mechanism for intermediate-depth earthquakes in subducting slabs and continental root zones.

Figure 1. Comparison of field observations with the numerical simulations. The photomicrographs show a typical shear zone (A1) and a pseudotachylyte (A2) from the Kråkenes gabbro. The shear zone shows no signs of melting, whereas the pseudotachylyte has a molten core and limited deformation at the margins. Rows B-E show the results of two numerical simulations that reproduce the observed structures. The left column illustrates a shear zone forming simulation and the right column a pseudotachylyte forming simulation. Both columns show the evolution of passive markers together with regions of complete melting (white) for model times shown in the lower plot, where the time dependence of the stress is illustrated. The predicted unloading of simulation 1 (D1) and the thermal relaxation in simulation 2 (E2) are fully consistent with the strain marker patterns recorded in the field. The white curves in A1 are initially vertical lines with the deformation field from D1.

References
Braeck, S. & Podladchikov, Y. Y., 2007 Spontaneous thermal runaway as an ultimate failure mechanism of materials. Phys. Rev. Letters 98.

John, T., S. Medvedev, L. Rüpke, T. B. Andersen Y. Podladchikov, and H. Austrheim, 2009. Generation of intermediate-depth earthquakes by self-localizing thermal runaway. Nature Geoscience, 2, 137-140.

Braeck, S., Y. Podladchikov, and S. Medvedev, (submitted). Spontaneous dissipation of elastic energy by self-localasing thermal runaway, Phys. Rev. E.

Top

Small-scale convection in the asthenosphere: A case study of foredeep basin formation
Elena Timoshkina and Valentin Mikhailov Institute of physics of the Earth, Russian academy of Sciences, Moscow, mikh@ifz.ru

We present detailed investigation of small-scale convection in the asthenosphere based of further developed thermo-mechanical model of evolution of rheologically stratified Earth outer shell. This shell includes four layers: sedimentary cover, lithosphere, asthenosphere and underlying upper mantle [1, 2]. The lithosphere - asthenosphere boundary is rheological one. Equation for the top of the model includes detailed description of sedimentation and erosion processes. The model permits modelling of active extension-compression by intraplate or mantle-induced forces as well as modelling of adjustment of mechanical and thermal equilibrium disturbed at active stages. The model demonstrates that after the Earth’s outer shell was disturbed by external intraplate or mantle-induced forces, small-scale convection within the asthenosphere arise. This convection makes for deformations in the lithosphere over a long period of time after cessation of external tectonic forces. In a compressional environment the small scale convection amplifies uplift of orogenic belts and causes subsidence at their periphery. We consider the small scale convection to be the main driving mechanisms of foredeep basins formation [3].

To illustrate this model we perform results of detailed modelling of the Great Caucasus orogen formation. To assign correctly initial conditions to the beginning of compressional stage, we considered preceding stages including: (1) extension of continental lithosphere in the early Jurassic; (2) subsequent post-extensional subsidence; (3) compressional (collisional) stage, when the system orogen – foredeeps forms. Parameters of the lithosphere and the asthenosphere and parameters of extensional – compressional processes were selected to provide a result close to the Great Caucasus - North Caucasus foredeeps, including topography, deep structure, thermal regime, subsidence history, gravity anomalies and so on.

The suggested model showed a good agreement with the data on the foredeeps structure and evolution. In particular, it is able to explain thickness of sediments in foredeep basins and their shape, formation of foredeeps not only at the front but also at the back of compressional thrust belts, uplift of a foredeep during compression in the belt and rapid subsidence after cessation of external compression.

Comparison of the numerical results with the observed data on the North Caucasus foredeep permits new interpretation of existing geological data [4]. In particular, it is possible to conclude, that the system orogen – foredeep was formed in a result of at least five active compressional events separated by periods of relatively weak tectonic activity. The first compressional event took place before the formation of the Maykopian series, i.e. 39.5 Ma, and could be related to the closure of the Arabian Ocean and subsequent beginning of the continent-continent collision in the Lesser Caucasus. There is still now consensus on when orogeny in the Caucasus region commenced, many researchers estimate beginning of the orogen formation by considerably later date. The four further compressional events can also be recognised - one of them being between 16.6 and 15.8 Ma, the others - between 14.3 and 13.7 Ma and between 7.0 and 5.2 Ma. These stages coincide well with geological data. The present day stage is also an active compression one.

Fig.1 Distribution of temperature (gray scale), position of main boundaries in the lithosphere and sedimentary cover and small-scale convection in the model of the Great Caucasus – North Caucasus foredeep formed in result of four compressional events before beginning of the present-day compressional one. The right side of the symmetric figure is shown. The centre of the orogen is at the left (x=0). The maximum length of arrows corresponds to 1.3 mm/year.

References
1. Mikhailov V.O., Myasnikov V.P., Timoshkina E.P. 1996 Dynamics of the Earth' outer shell evolution under extension and compression. (Izvestiya) Physics of the Solid Earth, v.32, N.6, pp. 496-502.
2. Mikhailov V.O., Myasnikov V.P., Timoshkina E.P. 1993 On the interaction of mantle flow with rheologically stratified boundary layer. Proceedings(Doklady) of Russian Academy of Sciences, v.330 (N6). p. 771-773.
3. Mikhailov V.O., Timoshkina E.P., Polino R. 1999 Foredeep basins: the main features and model of formation. Tectonophysics v. 308, pp. 345-360.
4. Mikhailov V. O., V. M. Gordin, E. P. Timoshkina, E. A. Kiseleva, and E. I. Smolyaninova 2007. Geodynamic Models and Their Application in the Combined Interpretation of Geological and Geophysical Data (a Review). (Izvestiya) Physics of the Solid Earth, v.43, N1, pp. 2-12.

Докладчик: Тимошкина Елена Павловна, к.ф-м-н., снс лаб. 501.

Top

Application of adaptive wavelets methods in computational geodynamics.
Yury Mishin (yury.mishin@erdw.ethz.ch)¹, Oleg Vasilyev² , Taras Gerya¹

1. Department of Earth Sciences, ETH Zurich, Switzerland
2. Department of Mechanical Engineering, University of Colorado, Boulder, USA

Calculation of phase diagrams for the Earth’s materials exerts a growing impact in studying of core and mantle dynamics. Recently, the method for calculation of phase diagrams and related in situ properties of rocks (e.g. density, enthalpy, seismic velocity, etc.) in space of two variables (P, T) by using second generation 2D wavelets was proposed (Vasilyev et al., EPSL, 2004). Here we present further generalization of this approach to three dimensions. We employ combination of adaptive wavelet-based meshing technology and efficient ‘‘phase diagram function’’ designed as a Gibbs free energy minimization. The proposed automated strategy allows one to obtain equilibrium phase assemblages and related physical properties depending on three arbitrary variables (e.g. P, T, water activity). The use of this strategy captures very small details of phase diagram morphology allowing both acceleration of the calculations and efficient compression of the results. Compressed phase diagrams can be efficiently used in geodynamic numerical modelling studies of tectonic processes on the basis of coupled petrological-thermomechanical modelling approach.

Top

3D numerical experiments assessing the necessary conditions for a plume-fed asthenosphere
Jason Phipps Morgan¹, Chao Shi¹, and Joerg Hasenclever²
1. EAS Dept., Cornell Univ., Ithaca, NY 14853 USA
2. Institute of Geophysics, Uni. Hamburg, Hamburg 20146 DE

In the past years we have presented observational evidence which suggests to us that in Earth’s mantle there exists a buoyant asthenosphere layer fed by upwelling in mantle plumes, and consumed by accretion and transformation into overlying lithosphere by ridge upwelling and melt-extraction (which creates a ~60km-thick layer of compositional lithosphere at mid-ocean ridges), by plate cooling (which accretes a further ~40km of asthenosphere after 100 Ma of near-surface cooling), and by dragdown by subducting slabs (which drags a further ~20km sheet of buoyant asthenosphere on either side of the subducting slab). This scenario has been recently reviewed in Yamamoto et al (GSA Vol. 431).

We believe that the reason this mode of mantle convection has not yet been seen in numerical models of mantle convection is due to the inability of current models to model the correct upwelling rates in focused lower-viscosity plumes (i.e. that, due to numerical resolution problems they currently underpredict plume upwelling) and to correctly model the magnitude of downdragging of a more buoyant but lower viscosity asthenosphere layer by subducting slabs (which they currently overpredict, cf. Phipps Morgan et al., Terra Nova, 2007). This happens because there are serious technical issues in accurately detemining flow with strong and correlated buoyancy and viscosity.

Here we present results from 3D calculations that include the effects of ridge accretion, plate cooling and well-resolved asthenosphere dragdown by subducting slabs. For our initial experiments, we place a ridge-centered plume at one corner of the computational region to take advantage of the two resulting symmetry planes that allow us to model the plume with only the nodes needed for similar resolution for a plume within the interior of the 3D computational region. We use an unstructured grid with high grid resolution within the plume conduit and slab entrainment regions of the box. We find that in order to create a global sub-oceanic plume-fed asthenosphere: 1) the plume-flux should be more than about 1.2 times the slab-flux; 2) the asthenosphere should be at least 0.5% more buoyant than underlying mantle when its viscosity is ~ 10e19 Pa-S, and ~0.05% more buoyant when its viscosity is 1e18 Pa-s.

We will also discuss in detail some of the numerical issues that arise during parallel iterative solutions for Stokes Flow with varying viscosity and buoyancy. We have found that the accuracy of the outer (Pressure Schur-complement) loop is contaminated in a curious way by the tolerance used to solve for the velocity subproblems, and have found a workaround to correct for this contamination. We also compare several iterative preconditioners, including the currently popular BfBt (also called LSC) approach.

Top

The Emergence of Low Lids in SuperSized Earths INFLUENCE OF THE PRESSURE ON THE CONVECTIVE BEHAVIOUR OF EARTHLIKE PLANETS
Lena Noack, Vlada Stamenković, and Doris Breuer Joint Planetary Interior Physics Research Group of the University Münster and Institute of Planetary Research DLR Berlin, Germany (lena.noack@dlr.de)

The convection pattern and the heat transport in a terrestrial planet depend on the viscosity of the mantle material, which defines the flowability of the material. The viscosity on the other hand depends strongly on temperature and pressure, i.e., the viscosity decreases with increasing temperature but increases with increasing pressure. Thus, the larger the planet the stronger is the influence of the pressure on the viscosity and a stagnant lid at the CMB can arise. To distinguish this layer from the lid at the surface, we use the notations lowlid and uplid [1].
In the present study, we used 2D and 3D spherical convection models as well as a parameterized model to examine the influence of the pressure on the viscosity, convection velocity, lid development, and heat flux at the CMB for Earthsize planets and SuperEarths of 2.5 Earth masses.

Simulation Model: The simulation runs are done using the 3Dspherical code GaiaS [2,3] as well as a 1D parameterized model based on thermal boundary layer theory [4]. The parameters of the Earthlike model [5] are rescaled to a SuperEarth planet of 2.5 Earth masses using the scaling laws of [6]

Viscosity Law: We use the Arrhenius law where the temperaturedependence is controlled using an activation energy E* and the pressureinfluence depends on the activation volume V*. In the model runs, we have chosen E* to be 150 kJ/mol and 300 kJ/mol.
Observational studies of glacial rebound on Earth suggest a maximum value for the activation volume of V*=2.5 cm³/mol [7]. However, experiments and theory suggest that the value might be much higher [8]. In the present study the activation volume ranges from V*=2.5 cm³/mol up to V*=5 cm³/mol.

Results: The higher the activation volume V*, the higher is the viscosity at the coremantle boundary. A high viscosity leads to an inner stagnant lid, the lowlid, and restricts the cooling of the core. Since the pressure at the CMB is larger for larger planets, the influence of the activation volume on the viscosity is much stronger for planets of 2.5 Earth masses than for Earth itself.

References
[1] Noack L, V. Stamenković, and D. Breuer (2009). ESLAB 09, P1.04.
[2] Hüttig C. and K. Stemmer (2008). Geochem. Geophys. Geosyst., 9 (2).
[3] Hüttig C. and K. Stemmer (2008). PEPI, 171, 137146.
[4] Grott, M. and D. Breuer (2008). Icarus, 193, 503 – 515.
[5] de Smet, J.H. (1999). Geologica Ultraiectina, No. 172.
[6] Valencia, D., R.J. O’Connell, and D.D. Sasselov (2006). Icarus, 181, 545554.
[7] Davies, G.F. (1999). Cambridge University Press, Cambridge.
[8] Karato, S. and P. Wu (1993). Science, 260, 771778.

Top

Novel methods in computational mineral physics: insights into minerals and compounds under pressure.
Artem R. Oganov
Dept. of Geosciences, Dept. of Physics and Astronomy, and New York Center for Computational Sciences, Stony Brook University, Stony Brook NY 11974-2100, USA.
Geology Department, Moscow State University, 119992 Moscow, Russia

In this talk I shall discuss some of the breakthroughs in computational mineral physics that occurred in the last several years. These advances put mineral physics on an altogether new level and include:

1. Ab initio solution of the crystal structure prediction problem using the USPEX method developed by us [1]. Crystal structure prediction has long been considered an insoluble problem, but this new method turns out to solve this problem very efficiently. I shall discuss some of the applications of this method, including prediction of new minerals [2] and novel phenomena under pressure – e.g. partially ionic form of boron [3] possessing superhardness [4] and a transparent phase of sodium under pressure [5].

2. First applications [6] of the ab initio metadynamics methodology to geophysically relevant problems that include crystal structure prediction, studies of phase transition and plastic deformation mechanisms of mantle minerals – MgSiO3 perovskite and post-perovskite.

3. The possibility to calculate mineral phase diagrams from first principles, as has been done by us for MgO, SiO2, Al2O3, MgSiO3 (e.g. [7-9]).

4. Discovery and intense studies of MgSiO3 post-perovskite [9,10], the main mineral of the D” layer.

5. Studies of the elastic constants of minerals (e.g. [11]) and their implications for the temperature distribution in the Earth.

1. Oganov A.R., Glass C.W. (2006). Crystal structure prediction using ab initio evolutionary techniques: principles and applications. J. Chem. Phys. 124, art. 244704.
2. Oganov A.R., et al. (2008). Novel high-pressure structures of MgCO3, CaCO3 and CO2 and their role in the Earth’s lower mantle. Earth Planet. Sci. Lett. 273, 38-47.
3. Oganov A.R., et al. (2009). Ionic high-pressure form of elemental boron. Nature 457, 863-867.
4. Solozhenko V.L., Kurakevych O.O., Oganov A.R. (2008). On the hardness of a new boron phase, orthorhombic γ-B28. J. Superhard Mater. 30, 428-429.
5. Ma Y., Eremets M.I., Oganov A.R., et al. (2009). Transparent dense sodium. Nature 458, 182-185.
6. Oganov A.R., Martoňák R., Laio A., Raiteri P., Parrinello M. (2005). Anisotropy of Earth’s D” layer and stacking faults in the MgSiO3 post-perovskite phase. Nature 438, 1142-1144.
7. Oganov A.R., Gillan M.J., Price G.D. (2005). Structural stability of silica at high pressures and temperatures. Phys. Rev. B71, art. 064104.
8. Oganov A.R., S. Ono (2005). The high-pressure phase of alumina and implications for Earth’s D” layer. Proc. Natl. Acad. Sci. 102, 10828-10831.
9. Oganov A.R., Ono S. (2004). Theoretical and experimental evidence for a post-perovskite phase of MgSiO3 in Earth's D” layer. Nature 430, 445-448.
10. Murakami M., Hirose K., Kawamura K., Sata N., Ohishi Y. (2004). Post-perovskite phase transition in MgSiO3. Science 304, 855-858.
11. Oganov A.R., Brodholt J.P., Price G.D. (2001). The elastic constants of MgSiO3 perovskite at pressures and temperatures of the Earth's mantle. Nature 411, 934-937.

Top

Collisional granitoids formation in rheologically layered lithosphere of overthrusted structures – numerical modelin
O.I. Parphenuk Institute of Physics of the Earth RAS, Russia, oparfenuk@mail.ru

Continental collisional structures are characterized by the crust thickening, high-grade metamorphic rocks and marked by positive gravity and aeromagnetic anomalies. These fundamental common features reflect the effect of the main tectonic event – horizontal shortening in compressional setting, collision of two continental plates by the mechanism of upthrusting accompanied by the thickening of the crust and the surface uplift. The main petrological mark of such collision is granite melts generated at different depth’s levels and exposed at the surface of the terranes and collisional shear zones.

Deeply eroded areas of the Archean and Proterozoic continental shields formed in the process of multistage tectonic evolution including horizontal shortening and collision expose at the surface middle to the lower crustal rocks uplifted by overthrusting from the depths of 20-40 km [1,2]. Extensive development of horizontal and oblique motions of crustal plates and blocks leads to nonstationary disturbances in the thermal regime, the heat flow and the surface and Moho topography.

Thermal-mechanical models of continental collision including horizontal shortening, brittle overthrusting in the upper crust and lower crustal viscous flow are applied to the simulation of the thermal and tectonic evolution of the overthrusted structures. Finite-element 2-D modeling was used to examine the conditions under which ductile flow of the rheologically layered lower crust and upper mantle can produce the gravitationally unstable structures with crustal roots and surface uplift as a result of shortening, nonuniform loading and erosion [3]. The thermal effects of collisional process and postorogenic stage are studied.

Thermal calculations show that the main temperature increase is observed at depths of the middle and lower crust and is fairly significant (up to 250K). In other words, the temperature characteristic of depths of 40-60 km is attained at depths of 20-40 km, creating the conditions for the partial melting and granite melt generation [4]. This observation confirms the idea of crustal origin of the continental collisional granites [5]. At the deeper lithospheric mantle levels the collisional geotherms mainly follow the deformation and are primarily controlled by initial conditions.

Since the model is predominantly two-dimensional, it can provide insights into the variety of metamorphic conditions in a region that experiences deformations in a horizontal compression setting. Collisional model confirms the observation of a wide range of P - T conditions over short distances in metamorphic belts: as a result of upward motion along the fault, the additional loading redistribution in the course of erosion, and the viscous compensation at the level of the lower crust P - T histories will be completely different for the points along the overthrust fault. For example, the material will be transported from a depth of 20 km to a depth of about 4.5 km while the rocks that were initially located at a depth of 3.5 km will experience a rather complicated P-T evolution [4]. The heat flow drops by about 10 mW/m2 above the bounding fault during obduction. The early postcollisional stages are characterized by the significant heat flow increase because of the thickening of the heat producing upper crust layer.

The very slow uplift and resulting erosion of the upper crust at collisional and postcollisional stages led to a reduced crustal heat production and to the increased strength of the lithosphere. These processes can result in the preservation of the crustal roots in such Precambrian structures as the Kapuskasing uplift (the Canadian shield) and collisional zones of the Anabar shield (the Siberian craton) for 2 Ga [1-3].

The sedimentary basins are formed at the edges of the collisional uplift by the rocks denudated from the surface elevation. At the frontal edge of the thrust sheet the basin depth reaches the value of more than 4 km with the maximum erosion level of ~9 km of uplifted crust. At the postcollisional stage the compressional regime is changed by extension but under the stresses of an order of magnitude less. The rate of the uplift depends on the viscosity values and the rate of erosion.

Significant local effect may result from the frictional heating along the slip zone during the overthrusting [6]. The frictional heat production is proportional to the product of the shear stress across the slip plane, the velocity of obduction and coefficient of friction. The additional heating can rise the temperature by 50-150ºC at 10-20 Ma in the vicinity of slip zone in the case of the rate of horizontal shortening up to 4 cm/yr and thrust sheet thickness up to 20 km. But this local and moderate additional heating (in comparison with crust thickening effect) may initiate partial melting process at depth levels of the middle crust [6,7]. The massive layer of granite melt approximately of 10 km thickness is observed at 10-15 km depth by geophysical methods in the resent collisional systems, such as Himalayas and Caucasus [1]. If this happens the produced frictional heat sharply lessen. As the relaxation of the local thermal anomalies proceeds rapidly enough, the melting may become a self-limiting process and overthrusting will take place with variable velocity.

The research is supported by the Russian Foundation for the Basic Research (grants 06-05-65221, 09-05-01032).

References
1. Rosen O.M., Fedorovsky V.S. Collisional granitoids and the Earth crust layering. Moscow, Nauchnyi Mir, 2001. 188 p. (in Russian).
2. Perry H.K.C., Mareschal J.-C., Jaupart C. Variations of strength and loсalized deformation in cratons: The 1.9 Ga Kapuskasing uplift, Superior Province, Canada // Earth Planet. Sci. Lett. 2006. V. 249. P. 216 – 228.
3. Parphenuk O.I., Dechoux V., Mareschal J.-C. Finite-element models of evolution for the Kapuskasing structural zone // Can. J. Earth Sci. 1994. V. 31. P. 1227 – 1234.
4. Parphenuk O.I. Thermal regime of collisional overthrust structures // Izvestiya, Physics of the Solid Earth. 2005. V. 41, No 3. P. 68 – 70.
5. Hain V.E., Lobkovskii L.I. Some features of the collisional orogens formation // Geotectonics. 1990. V. 6. P. 20 – 31 (in Russian).
6. Brewer J. Thermal effects of thrust faulting // Earth Planet. Sci. Lett. 1981. V. 56. P. 233 – 244.
7. Cermak V., Bodri L. Time-dependent crustal temperature modeling: Central Alps // Tectonophysics. 1996. V. 257. P. 7 – 24.

Top

Modeling Transform Plate Boundary in 3D: Dead Sea Transform Fault from the Red Sea to Lebanon Mountains
Alexey G. Petrunin^1,2, Stephan V. Sobolev1 and Zvi Garfunke^1,3
1. GeoForschungsZentrum, Potsdam, Germany (alexei@gfz-potsdam.de)
2. Institute of Physics of the Earth, Moscow, Russia
3. Institute of Earth Sciences, Hebrew University of Jerusalem, Israel (zvi.garfunkel@huji.ac.il)

Dead Sea Transform (DST) fault system is a part of the Syrian-African rift system and it extends from the divergent plate boundary of Red Sea rift at the south to the convergent plate boundary in the Taurus Mountains at the north. The DST itself is a left-lateral transform fault, accommodating differential motion between the African and Arabian plates. The morphology of the DST fault system is expressed by several linear stretches separated by a number of pull-apart basins, where the Dead Sea basin is a largest. Our previous models (Sobolev et al. 2005, Petrunin and Sobolev, 2006, 2008) have been focused at two main topics: (1) major controls of the fault localization in strike slip settings and (2) major controls of the structure and evolution of pull-apart basins located at strike- slip faults. To do so, we modeled lithospheric deformations using 2D and 3D FEM techniques with realistic elasto-visco-plastic temperature and stress dependent rheology. However, limitation of our models was their relatively small size, which did not allow including the source of the strike-slip motion in the region (opening of the Red Sea Rift), and major obstacle for the propagating fault resulting in its bending in Lebanon. In present work we extend the model to the larger region. The new model domain includes north-western part of the Red Sea and extends to the Lebanon Mountains in the north where deformation becomes more complicated and large part of the strike-slip motion becomes distributed. Because of the large size of the domain, we made an improvement in the modeling technique taking into account sphericity of the Earth surface.

Here we present a model aimed at revealing controls of initiation and positioning of the DST fault. We also consider the discrepancy between present day surface heat flow (about 50mW/m2, consistent with the thickness of the lithosphere of more than 120 km) and observed thickness of the lithosphere at the DST area of 70-80 km according to seismic data (Mohsen et al. 2005). This discrepancy means that lithosphere around DST was thinned in the past and related high heat flow had not enough time to reach the surface.

As an initial setup we use simplified geometry of the lithosphere-asthenosphere boundary and lithospheric structure corresponding to continental margin conditions. The lithospheric thickness is reduced by 50 km at some time, which is considered as a model parameter to fit observations. In contrast to previous models, where kinematic boundary conditions were used (constant velocities at the side boundaries), in this model we apply a constant force at a right boundary of the modeling domain, which is also considered as model parameter. The model shows evolution of strain distribution during the Red Sea opening. As the Red Sea basin spreads, the strain is first accumulated in broad zone along the DST during several millions years. The strike-slip rate and displacements at this stage are low, less than 20% of their present-day values. After thinning of the lithosphere by about 50 km the strike-slip strain localizes at the region of minimum strength controlled by crustal thickness and thermal structure of the lithosphere.

From our model we conclude that to localize the strike slip deformation at DST in the cold lithosphere (surface heat flow lower than 50mW/m2) the following conditions had to be met: (i) large force (about 1.6x1013N/m) applied by the Arabian plate, (ii) relatively low strength of mantle lithosphere, (iii) thermal erosion of the lithosphere that triggered the DST localization.

Top

Pressure variations during ultra-high pressure metasomatism
Y. Y. Podladchikov and J. C. Vrijmoed PGP, University of Oslo, Physics of Geological Processes (PGP), Oslo, Norway (yuripo@fys.uio.no)

Incomplete reactions or apparent disequilibrium between minerals in metamorphic rocks from high pressure (HP) and ultra-high pressure (UHP) terrains are commonly observed in field and thin section. Here we study an example of a peridotitic body enclosed in migmatitic felsic gneiss at the UHP terrain of the Western Gneiss Region (WGR) in Western Norway. The WGR represents the basement of Baltica which became metamorphosed during the Caledonian Orogeny. In the highest grade parts of the WGR metamorphic temperatures were between 600-800 C for several Ma, pressures reached the diamond stability field (5 GPa), fluids were available and yet well known pressure sensitive reactions such as from spinel- to garnet-peridotite went only locally to completion in the studied peridotite. Previously reported field observations, major, trace and isotope geochemistry, geochronology, mineral-chemistry and textures suggest that this peridotite body became metasomatised by supercritical fluids from the host rock gneiss during the Caledonian Orogeny. Contrasting pressure estimates, incomplete reactions and preserved compositional gradients in minerals from this body may indicate very rapid exhumation, very sluggish kinetics of diffusion, metasomatically disturbed equilibrium at low temperatures or a combination of the three. These possibilities call for a re-evaluation of existing methods in geothermobarometry. Alternatively we may consider that pressure variations existed from the grain up to the outcrop scale during the metasomatism of the peridotite at UHP conditions. The mechanical responses to volume changes that are involved in chemical reactions in rocks may control the progress of reactions significantly (Schmid et al.1, 2009). Our main results, obtained from numerical modelling, show that pressure variations may be generated, and maintained on the geological time scale, as a result of the mechanical response of the rocks during melting reactions in confined space. Geodynamics implications of our alternative explanation will be discussed.

Top

Numerical modeling of Rayleigh-Tailor instability in relation to the melting and diapirism of granite magma in the crust
Polyansky O. P.¹, Korobeynikov S. N.², Babichev A. V.¹, Reverdatto V. V.¹
1. Institute of Geology and Mineralogy, Novosibirsk, pol@uiggm.nsc.ru
2. Lavrentyev Institute of Hydrodynamics, Novosibirsk, korob@hydro.nsc.ru

Method for mathematical modeling of partial melting and ensuing development of gravitational instability of acid melt in the thickened earth crust with enlarged thickness of granite layer has been suggested.

Statement of the problem is as follows:
1) 2D right-angled area in the earth crust measured 38x60 km;

2) The equations for mechanical equilibrium in “weak” form were solved (equation of the principle of virtual work);

3) The heat transfer equation was solved with consideration for phase transition on melting and with the use of variable coefficients;

4) Instead of classical statement of the Stefan problem the equation having regard to the increased heat capacity at phase transition was put to use;

5) The governing equations of quasi-static deformation were solved numerically within the framework of an approximation of plane stain problem;

6) The finite elements method was used for the numerical discretization of equations of deformable solid continuum mechanics (DSCM);

7) The numerical modeling was performed using the commercial MSC.Marc software in which all types of non-linearity of the DSCM were taken into account.
The model has been constructed for the investigation of melting and the buoyant upwelling of granitic magma caused by this process. The heating was induced by underplating of basic magma in the base of continental crust. The purpose of modeling is to estimate the parameters of melting and diapirism of acid melt in the lower crust, to determinate the flow structure of buoyant granite magma and to predict possible form of granite-gneiss bodies.
The first set of experiments has been carried out at constant yield stress σY= 1 MPa [1] which is temperature-independent. The results of modeling are displayed in Fig.1a-c in the form of temperature field for the final state of melt upwelling. At the beginning, the propagation of partial melting above the heat source, i.e. basic intrusion, took place in the course of c. 2 Myr. The conductive heating which proceeded in the absence of motion of crustal material has gone on throughout c. 2 Myr at given parameters. In this period the melted zone was developed in dome shape 6-7 km in height and at given width of basic intrusive body (20 km). Thereafter the upwelling of melted material has begun (Fig. 1a; total time span is 2.162 Myr).
The second set of calculation experiments has been performed for revealing the relationship between the buoyant upwelling dynamics of diapir and the temperature-dependence of plasticity. Fig. 1b presents the results of calculations for temperature-dependent plasticity when the yield stress is variable from 10 to1 MPa in the range from 20°C (at the surface) to 650 °C (melting). Compared to the previous set of calculations, the ascent of diapir was bounded by the level of middle crust at a depth of 14-16 km (Fig. 1b). Time span is 7.95 Myr.
Fig. 1c presents the results of modelling for “softer” rheology of the crust: yeld stress is variable from 5 to 0.5 MPa. Two branches of upwellig diapir are developed. Total time span is 2.132 Myr.

The numerical experiments allow to draw the conclusions:
1) Critical volume of partial molten mass is bound to appear before the buoyant upwelling can begin in gravitational field. According to model estimates, the height of molten area within the granite crust is due to be no less than 6-7 km.

2) Independently of size of the heat source (with fixed or variable width), the mushroom-shaped upwelling bodies occur in all model versions. The incurrent canal of high-temperature melt and the head of diapir are distinguished.

3) The height of diapir upwelling depends on the rheology of surrounding rocks. An increase of yield stress limit of the order (from 1 to 10 MPa) at decreasing temperature restricts a feasible level of uplifting to a depth of 10-15 km. The similar estimates have been obtained when modeling diapirism with the supposition of viscous rheological behavior of lower-crustal material [2]. Above the axial part of diapir the crustal surface uplift forms. A geological example of the model is the Teya granite-gneiss dome located in the Trans-Angarian region at southwestern margin of Siberia platform. Its area is more than 2000 km² [3]. It occurs in axial part of the Central Yenisey ridge uplift, made up of metamorphic and granitized rocks. Large-scale geomapping and petrogeochemical study indicate that the dome core is made up mainly of gneisses, porphyroblastic gneisses, granite-gneisses and intrusive granites in mean density of 2600 kg/m3. The country rocks are made up of biotite-bearing phyllites, quartz-mica ±Gr±St±And±Sil schists in density of 2960 kg/m3 and orthoamfibolites in density of 2920 kg/m³.
The determination of U-Pb age with zircon indicates that one of emplacement stage of the Teya granite-gneiss dome has occurred at 866±16 Myr. More recently, this age was confirmed by SHRIMP analysis of zircon grains, which were sampled from granitoid rocks in northern part of the Teya massif. Our model calculations have pointed to the faster ascent velocity of melt body in relation to rate of the conductive heating and cooling front of granite massif. According to model calculation the duration of buoyant upwelling, cooling and solidification of granite diapirs covers no more than 7-8 Myr. This estimate is less than the accuracy of isotope geochronological determinations. Hence the model estimates allow to interpret the obtained geochronological data as an age of ascent and solidification of granite melt “at the site”.

1. Gerya T. V., Burg J.-P. Intrusion of ultramafic magmatic bodies into the continental crust: numerical solution // Phys. Earth Planet. Interiors. 2007. V.160. P.124-142.

2. Bittner D., Schmeling H. Numerical modelling of melting processes and induced diapirism in the lower crust // Geophys. J. Int. 1995. V. 123. P.59-70.

3. Nozhkin A. D., Turkina O. M., Bibikova E. V., Terleev A. A., and Khomentovsky V. V. Riphean granite-gneiss domes of the Yenisei range: geological structure and U-Pb isotopic age// Russian Geology and Geophysics. 1999. Vol. 40, No. 9, pp. 1284-1292.

Fig1. The results of modeling are displayed in the form of temperature fields for the final states of melt upwelling under constant (a) and temperature-dependent yield stress (b-c). See comments in the text.

Top

StGermain software environment for HPC and software maintenance.
S Quenette, L Hodkinson, L Moresi, B Appelbe

Around the world there are initiatives to consolidate the effort expended on creating infrastructure for science. For example with respect to modelling software infrastructure for geophysics, in Australia there is an initiative named AuScope Simulation & Modelling, in the USA the equivalent initiative is CIG, and so on.

The fundamental challenges are to combine the utilisation of ever changing High Performance Computing, with the ever-increasing complexity in numerics and physics, whilst appropriately dealing with the subsequent explosion in software complexity.

Figure 1 - This is an example of AuScope supporting research infrastructure. The image is a test case for a basin evolution model in Underworld coupled to a surface process model in SPModel. Both codes are HPC distributed memory parallel, and extensible to many applications.

This scenario is a catch-22 for geophysics researchers. On one hand, this increasing "complexity" is rendering the chance of a home-grown code leading to disruptive science as increasingly less likely. On the other hand, a community code may have a researcher's contribution "lost within it", "impractical to merge with", or it may be too difficult to learn.

Our software environment - StGermain, and the derivative Underworld framework and examples like the GALE application, aims to address these problems. In doing so, it provides an environment for the contribution of scientific concerns (whether they be numerical, physics based or geological phenomena) by a community of multidisciplinary researchers, in a software maintainable manner. Discussed are the approach, what we've achieved thus far, and where we are headed.

Top

Magnetic filed generation in Earth’s core
Reshetnyak M. Institute of the Physics of the Earth RAS, Moscow, e:mail: m.reshetnyak@gmail.com

Many astrophysical objects such as galaxies, stars, the Earth, and some planets have large-scale magnetic fields that are believed to be generated by a common universal mechanism - the conversion of kinetic energy into magnetic energy in a turbulent rotating shell. The details, however, and thus the nature of the resulting field, differ greatly. The challenge for the dynamo theory, see, e.g., [1], is to provide a model that can explain the visible features of the field with realistic assumptions of the model parameters. Due to the strong non-linearity of the planetary dynamo regimes the usual approach is a 3D direct numerical simulations (DNS) which try to mimic processes of the thermal and compositional convection in the liquid metal core and the self-generation of the magnetic field in the conductive medium. Calculations for the entire planet are done using either spectral models or finite-volume methods and finite differences and have demonstrated beyond reasonable doubt that the turbulent 3D convection of the conductive fluid can generate a large-scale magnetic field similar to the one associated with small random fluctuations, see review in [2].

However, after some time when everybody excepted that even 3D models of geodynamo based on Boussinesq convection can reproduce spatial spectrum of the geomagnetic field, reversals and excursions, inner core rotation comparable with any desired observations, the question arised: what is the value of all these efforts, when the more or less the same morphology of the magnetic field at the planetary surface is produced by the very different mechanisms inside of the Earth? In other words the inverse geodynamo modelling (when the aim is to improve our knowledge on the invisible fields in the core and parameters in the models used) is an example of the very ill-poised problem. For the first 3D geodynamo models this question was really crucial, because without any assumptions and thorough justifications all diffusion coefficients were chosen to provide numerical stability of the model (e.g., the so-called hyperdiffusion) using some uncertain turbulent values. That is why some time required to analyse results of simulations and provide thorough research of the properties of the dynamo models in the spherical shells as a function on the diffusion coefficients, rotation rate of the shell and nature and magnitude of the driven forces [3,4].

Any way, the rigorous models of the turbulence for the liquid core was and is a challenge for the modern dynamo theory. So, even for the geodynamo (which is quite a modest case on the astrophysical scale) the hydrodynamic Reynolds number estimated on the west-drift velocity is Re ≈ 109, that requires N~Re^9/4 ≈ 10²⁰ modes to be resolved in the model, what is out of reach for any modern supercomputer. On the other hand the typical magnetic Reynolds number in the liquid core is in the range of 10²÷10³, what is already acceptable for the modern DNS, the best of which reach the mesh (or number of basic functions) up to 1000³. Even this resolution is in some sense is not necessary for geodynamo modelling because the accuracy of the measured by the satellites potential part of the geomagnetic field starting from 14^th spherical harmonics at the surface of the Earth falls to the first harmonics in paleomagnetic records already in some millions in the past. We come to dramatic situation: convection requires fine resolution, which makes simulations too slow if possible, and the long term paleomagnetic observations have very pure spatial resolution.

Obviously, the remedy is to develop suitable parameterisation of the small-scale fields what is a subject of the magnetohydrodynamic turbulence models. Whereupon the simplified large-scaled equations for the mean field should be solved. This technique is a subject of the mean field dynamo theory [5] was and is still used in many objects, e.g. Galaxy, Sun. The direct application of this approach to the Earth is complicated because of the rapid daily rotation of the planet (the Rossby number for the Earth’s core is Ro ≈ 10^-5). This geostrophic regime, which leads to the intermediate state inherited properties of 2D and 3D systems, can not be described by the commonly used Kolmogorov's like models of the turbulence. The level of anisotropy in such regimes can reach some orders of magnitude. The other point is possibility of the synenergetics effects even in pure convection for the rapid rotation. This effect for years was prescribed to the property of the averaged induction equation, which could describe generation of the large scale fields by the small-scale turbulence when the mirror-reflexional symmetry is broken (α-effect). At the same time from seventies was known that 2D Navier-Stokes equation itself provides inverse cascade of the kinetic energy (energy comes from the small scales to the mean scale) well-known in the physics of the atmosphere and ocean. It appears that in quasi-geostrophic turbulence inverse cascade of the kinetic energy holds at the scales where the magnetic field is generated [6]. Such an analysis of energy redistribution in k-space can be helpful for understanding of the nature of the small-scale physical processes covered by the mantle as well as for adjustment of the semi-empirical models of the turbulence for the realistic geophysical parameters. The anisotropy of the fields exists not only in the physical space, but comes to the wave space as well. As a result in contrast to the local energy transfer and local wave interactions predicted by the Kolmogorov theory of the isotropic and homogeneous turbulence, geostrophy leads to non-locality in the wave space. It means that the small scale motions can transfer energy directly to the large scales. This stands new requirements to the turbulent modelling in future.

References
1. Hollerbach, R., Rüdiger, R. The Magnetic Universe, Wiley-VCH Verlag GmbH & Co.KGaA, Weinheim, 2004. 338p.
2. Olson, P. ed. Core dynamics. Vol. 8 of Treatise on geophysics (Schubert, G. main ed.). Amsterdam: Elsevier. 2007. 358p.
3. Christensen, U.R., Olson, P., Glatzmaier, G.A. Numerical modeling of the geodynamo: A systematic parameter study // Geophys. J. Int. 1999. 138. P.393-409.
4. Christensen, U.R., Aubert, J. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields // Geophys. J. Int. 2006. V.166. P.97-114.
5. Krause, F., Rädler, K.-H. Mean field magnetohydrodynamics and dynamo theory. Berlin: Akademie-Verlag. 1980. 272p.
6. Reshetnyak, M., Hejda, P. Direct and inverse cascades in the geodynamo // Nonlin. Proc. Geophys. 2008. V.15. P.873-880.

Top

Modeling of stresses caused by lithospheric density inhomogeneities: delamination of Sierra-Nevada deep batholithic root
Romanyuk T.V. (Institute Physics of the Earth, Moscow, t.romanyuk@mail.ru)

In Southern California, beneath southern parts of Great Valley and Sierra-Nevada, a high seismic velocity body in the upper mantle is mapped by seismic tomography method [1, 2, 6]. It is called “mantle drip” and is similar to a cilindr with radia of ~60 km located between depths of ~50-250 км. The “mantle drip” is interpreted as a delamination of deep batholitic root of Sierra-Nevada batholith [3 -5, 7-11]. A 2-D density and rheological model for a profile “SN” cut the “mantle drip” over its center is compiled based on geological and geophysical data. 2-D calculations of stresses caused by density inhomogeneities are carried out for “SN” profile. A crusial moment for stress partern into mantle is a relation of rheology of the “mantle drip” and surounding mantle. If “mantle drip” is “harder” than mantle (Model «А»), then tangential stresses in the mantle are maximal beneath “mantle drip”. If “mantle drip” is “softer” than mantle (Model «B»), then tangential stresses in the mantle are maximal over the top and near the sides of “mantle drip”. Indirect arguments allow to conclude that Model “B” is preffered. A layer with overlithostatic extension is obtained in the middle-lower crust for Model “B”. The top of the layer is conside with lower limit of seismicity in the region. This layer plays a special role in migration of San Andreas fault system inside North American continent.

1. Biasi G.P., Humphreys E.D. P-wave image of the upper mantle structure of Central California and southern Nevada // Geoph. Res. Lett. 1992 V.19. p.1161–1165.
2. Boyd O.S., Jones C.H., Sheehan A.F. Foundering lithosphere imaged beneath the southern Sierra Nevada, California, USA // Science. 2004. V. 305, p. 660–662.
3. Ducea M.N., Saleeby J.B. A case for delamination of the deep batholithic crust beneath the Sierra Nevada // International Geology Review. 1998a. V. 40, p. 78-93.
4. Elkins-Tanton L.T., Grove T.L. 2003. Evidence for deep melting of hydrous metasomatized mantle: Pliocene highpotassium magmas from the Sierra Nevadas, J. Geophys. Res. V. 108 N B7, 2350, doi:10.1029/2002JB002168, 2003.
5. Farmer G.L., Glazner A.F., Manley C.R. Did lithospheric delamination trigger late Cenozoic potassic volcanism in the southern Sierra Nevada, California? // Geol. Soc. Am. Bull. 2002. V. 114, p. 754-768.
6. Fliedner M.M., Klemperer S.L., Christensen N.I. Three-dimensional seismic model of the Sierra Nevada arc, California, and its implications for crustal and upper mantle composition // J. Geoph. Res. 2000. V. 105, N. B5, p. 10899–10921.
7. Jones C.H., Farmer G. L., Unruh J. Tectonics of Pliocene delamination of lithosphere of the Sierra Nevada, California // Geol. Soc. Am. Bull. 2004. V.116. N.11-12. p.1408-1422.
8. Lee C.-T., R. L. Rudnick, G. Brimhall, Deep lithospheric dynamics beneath the Sierra Nevada during the Mesozoic and Cenozoic as inferred from xenolithpetrology , Geochem. Geophys. Geosyst., 2, 10.1029/ 2001GC000152, 2001.
9. Lee C.-T.A, Cheng X.., Horodyskyj U. 2006. The development and refinement of continental arcs by primary basaltic magmatism, garnet pyroxenite accumulation, basaltic recharge and delamination: insights from the Sierra Nevada // Contrib. Mineral. Petrol. V. 151. P.222-242. DOI 10.1007/s00410-005-0056-1.
10. Manley C.R., Glazner A.F., Farmer G.L. Timing of volcanism in the Sierra Nevada of California: Evidence for Pliocene delamination of the batholithic root? // Geology. 2000. V. 28, p. 811–814.
11. Zandt G., Gilbert H., Owens T., Ducea M., Saleeby J., Jones C. Active foundering of a continental arc root beneath the southern Sierra Nevada in California // Nature. 2004. V 431. p. 41-46.

Top

ROLE OF SEEPAGE FORCES ON FAULTING AND EARTHQUAKE TRIGGERING
Alexander Rozhko Weatherford Laboratories, Trondheimy, Norway

The earthquake triggering occurs as a result of the redistribution of stress induced by an earthquake or fluid pressure, while seepage forces caused by gradient in pore-fluid pressure [1]. In the recent publication [2] it was demonstrated how the location of aftershocks correlates with propagation of pore-fluid diffusion front. In this work we calculate the stress transfer caused by diffusion of pore-fluid using the analytical solution, derived for this purpose and demonstrate how the observed locations of aftershocks can be predicted, using the analytical solution for seepage forces.

Reference:
[1] Rozhko, A. Y., Y. Y. Podladchikov, and F. Renard (2007), Failure patterns caused by localized rise in pore-fluid overpressure and effective strength of rocks, Geophys. Res. Lett., 34, L22304, doi:10.1029/2007GL031696.
[2] Miller, S. A., C. Collettini, L. Chiaraluce, M. Cocco, M. Barchi, and B. Kaus (2004), Aftershocks driven by a high pressure CO2 source at depth, Nature, 427, 724– 727.

Top

Top

Calculation of site-specific carbon-isotope fractionation in pedogenic oxide minerals
James R. Rustad¹, Piotr Zarzycki¹, and Jean-Francois Boily²
1. University of California, Davis, CA USA
2. University of Umea, Sweden

The paleosol carbon-isotope method is one of the major techniques for reconstructing the composition of the ancient atmosphere. The depth-dependent mixing of respirative CO2 from soil (13C-depleted) with atmospheric CO2 (13C-enriched) can be preserved in soil minerals and has been used to infer atmospheric PCO2 and respiration rates in paleosols formed back into Precambrian times. The fractionation between mineral-dissolved CO2 [CO2(m)] and CO2(g) is a key factor in the technique because it determines how the respired CO2 vs. atmospheric CO2 mixing profile is imprinted onto the soil mineralogy. Ab initio molecular dynamics and quantum chemistry techniques are used to calculate the structure, vibrational frequencies, and carbon-isotope fractionation factors of the carbon dioxide component [CO2(m)] of soil (oxy)hydroxide minerals goethite, diaspore, and gibbsite. We have identified two possible pathways of incorporation of CO2(m) into (oxy)hydroxide crystal structures: one in which the C4+ substitutes for four H+ [CO2(m)A] and another in which C4+ substitutes for (Al3+,Fe3+) + H+ [CO2(m)B]. Calculations of isotope fractionation factors give large differences between the two structures, with the CO2(m)A being isotopically lighter than CO2(m)B by ≈10 per mil in the case of gibbsite and nearly 20 per mil in the case of goethite. boThe surprisingly large difference in the carbon-isotope fractionation factor between the CO2(m)A and CO2(m)B structures within a given mineral suggests that the isotopic signatures of soil (oxy)hydroxide could be heterogeneous. The calculated site heterogeneity is in good agreement with high-resolution 2-D correlation IR spectroscopy. The ability to probe site-specific isotope fractionation factors in minerals has the possibility to usher in a new era of high-resolution paleoclimatology.

Top

On the mantle control of core convection and the geomagnetic field
Ataru Sakuraba Department of Earth and Planetary Science, University of Tokyo, Tokyo 113-0033, Japan

Convection of the Earth’s fluid core and the resulting geomagnetic field are intrinsically controlled by the physical state of the mantle. Without a significant amount of radiogenic heat sources, thermal convection is impossible in the core if the mantle does not draw heat from inside. The mantle also controls compositional convection due to solidification of the solid inner core because the inner core growth needs overall cooling. The mantle control of the geodynamo has been mainly studied by imposing a laterally non-uniform heat flow at the core-mantle boundary in numerical modeling. It has been argued that the heat flow pattern affects frequency of geomagnetic reversals (Glatzmaier, Coe, Hongre and Roberts, 1999), geomagnetic secular variations like westward drift (Christensen and Olson, 2003) and the inner core growth (Aubert, Amit, Hulot and Olson, 2008). The heat flow pattern is often presumed by seismic velocity anomaly in the deepest part of the mantle. However, one-to-one correspondence between the seismic velocity and the temperature is highly arguable because of complicated compositional effects and phase transitions. The dynamical and thermochemical state of the deepest part of the mantle is one of the biggest issues in the geodynamo problem.

Because of computational reasons, previous numerical models of the geodynamo have imposed unrealistically large viscosity. Recent models have succeeded in reducing viscosity significantly. For example, Kageyama, Miyagoshi and Sato (2008) showed that reduction of viscosity led to small-scale, sheet-like convection by which a different dynamo mechanism ensued. Surprisingly, their solution exhibited a magnetic field not strongly dominated by an axial dipole field. This means that the solution becomes geophysically unrealistic though the parameter values approach an Earth-like regime. Paul Roberts and I recently demonstrated that this discrepancy arouse from improper use of thermal boundary condition (in submission). When the core surface temperature is laterally uniform, which is the same boundary condition as the Kageyama model, the convection pattern becomes small-scale and the generated magnetic field is relatively weak. The surface magnetic field is almost stationary and the field structure is broad. However, when the surface heat flux is forced to be laterally uniform, the convection pattern becomes large-scale and the magnetic field is dominated by a strong axial dipole (see Figure). There are several magnetic flux patches at the core surface and they move westward like the geomagnetic westward drift. We conclude that the isothermal boundary condition is not only geophysically unrealistic but also leads to failure to reproduce Earth-like strong magnetic fields even when more realistic parameter values are used.

I plan to investigate effects of laterally non-uniform heat flow in lower-viscosity geodynamo models. With a lower viscosity, we can reproduce magnetic field behaviors of small length scales and short time scales, so the numerical solutions can be compared to recent geomagnetic observations in detail. One of the robust geomagnetic field variations is the westward drift. It is well constrained by recent observational data that the westward drift is mainly seen in an equatorial belt in the Atlantic hemisphere and that the field variations are less prominent in the Pacific hemisphere. One possibility is that this hemispherical non-uniformity is caused by non-uniform cooling from the mantle. Taking account of ambiguity in interpreting the seismic anomaly near the bottom of the mantle by the temperature alone, I surmise that the geodynamo modeling could make a constraint on the thermal structure near the base of the mantle. That is, if a heat flow pattern is specified that explains best the geomagnetic secular variations, we could get information about the thermal structure independent of other seismic and geodetic data. I will present some results of this ongoing study and discuss how the geodynamo modeling can contribute to the study of the deep mantle.

Figure: Comparison of two different boundary conditions on the core surface temperature in a low-viscosity geodynamo model. The radial velocity on z = 0.1c is shown in each panel, where z is the distance from the equatorial plane and c is the core radius. The solution produced by an isothermal model is shown in the left, and the solution produced by a uniform heat flux model is shown in the right. All the parameters and the boundary conditions other then the thermal boundary condition at the core surface are the same in these two models. The velocity represents the magnetic Reynolds number.

Top

Density models and structure of the Earth`s lithosphere beneath Central and South Asia constrained with variations of free mantle surface
V.N. Senachin¹ and A.A. Baranov²
1. Institute of Marine Geology and Geophysics, FEB RAS, Yuzhno-Sakhalinsk, Russian Federation (geodyn@imgg.ru),
2. Schmidt Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russian Federation (baranov@ifz.ru)

Free mantle surface (FMS) is one of the important characteristics of the isostatic state of the Earth. FMS shows the degree of uplifting of the crust about the normal level, which corresponds to the homogeneous upper mantle. The FMS anomaly study can provide important information about the different geodynamic processes that responsible for the density heterogeneities in the upper mantle and changing isostatic state of the lithosphere. Investigations of the FMS (Artemjev et. al, 1983) revealed main dependencies for the depth of the FMS under the continents and oceans. For the continental lithosphere it was found that the FMS depth depends on the thickness of the crust. Subsequently, the same dependence was revealed for the oceanic lithosphere using CRUST 2.0 model for all Earth (Senachin, 2008). In this study we present the updated FMS anomaly map for the Central and Southern Asia (Fig.1) calculated using the crustal model AsCRUST-08 (Baranov, 2008), which has the resolution of 1x1 degree. We used the Moho map and density for upper, middle, and lower layers of crystalline crust for calculating the FMS anomalies.

Fig 1. Free mantle surface for Central and Southern Asia. Black border shows the border of As-Crust model, outside this area for calculations it was used Crust 2.0 model (Bassin et al, 2000)

The Southern and Central Asia is tectonically complex region characterized by the great collision between the Asian and Indian plates, anomalously thick uplifted crust, and the large extensional zones near the southern and eastern margins of Asia. The evolution of the entire region is also strongly related to the active subduction along the Pacific border. The crustal model AsCRUST-08 provides substantially more detailed FMS data for the Asia region. For example we can see anomalous uplifting of the FMS up to 3 km in the extensional zones (Red Sea) and in the deep seafloor areas. Arabian Peninsula has the FMS depth about 6 km, which can be attributed to rather high density of the upper mantle.

We have investigated density distribution of the upper mantle using FMS data. The main problem of modeling was to describe known connection between FMS and thickness of the Earth`s crust. Three models of lithosphere mantle have considered for define the reason of dependence: homogeneous mantle with constant density (model A), radial density variations in the mantle (model B), and lateral density variations in the mantle (model C).

Results of modeling have shown, that models A and B cannot explain observable FMS dependence on thickness of the crust; only model C gives quite suitable values of density. It is assumes the presence of high dense areas in isostatic compensated lithosphere which anomaly high density increases, with thickening of the crust.

References:
[1] Artemjev M. E., Chesnokov E M. Density inhomogeneities in the mantle using the free mantle surface data. Izvestiya, Physics of the Solid Earth, 1983, Vol. 5, p. 3-11.
[2] Baranov A. A. Complex 3D crustal model of Central and Southern Asia //Electronic journal «Vestnik RAS» №1(26)2008 (URL: http://www.scgis.ru/russian/cp1251/h_dgggms/1-2008/informbul-3_2008/cw-3.pdf (in Russian)
[3] Bassin C. Laske G., Masters G. The Current Limits of Resolution for Surface Wave Tomography in North America //EOS Trans AGU, 2000. - 81(48), Fall Meet. Suppl., Abstract F897. (http://mahi.ucsd.edu/Gabi/rem.html)
[4] Senachin V.N. The isostasy and density inhomogeneous of the lithosphere on the data of Crust 2.0 model //Relationship between the surface and deep structures of the earth’s crust: Proceedings of the 14th International Conference. Petrozavodsk: Karelian Research Centre, RAS. Part. 2. - 2008. - P. 186-188. (in Russian).

Top

Studying the Earth's interior based on correlations of ambient seismic noise
Nikolai Shapiro, Institut de Physique du Globe de Paris

Recent years witnessed a strong interest in studies based on background seismic noise. A large amount of continuously recorded digital seismic recordings provided by modern networks was used to develop new approaches for seismic tomography (e.g., Shapiro et al., 2005; Sabra et al. 2005; Yang et al., 2007) and monitoring (e.g., Sens-Schönfelder and Wegler, 2006; Brenguier et al., 2008a,b). Detailed analysis of high-quality continuous records also allows us to better understand the origin of the ambient seismic noise and its relation to the process in the ocean and atmosphere.

New methods of noise-based seismic imaging and monitoring are based on a principle that Green function between to points can be extracted by correlating a random wavefield recorded by receivers located in these points. In other words, one of two receivers can be considered as a virtual source recorded by the second receiver. This principle is especially attractive when applied in context of random wavefield recorded by a network of numerous recorders. In this case, by computing all possible inter-station cross-correlations it becomes eventually possible to place virtual sources at every receiver location and to have their records by the whole network resulting in very dense path coverage. The deterministic waveforms (Green functions) extracted from the cross-correlations can be used then to measure and to invert travel times with different methods developed in context of earthquake or explosion based seismology.

This possibility to extract deterministic Green functions from correlations of a random wavefield can be proved mathematically with different approaches (e.g., Snieder et al., 2004; Guedard et al., 2008) and has been demonstrated in acoustic laboratory experiments (e.g., Derode et al., 2003). It has been recently applied to extract surface-wave Green function from correlations of ambient seismic noise (e.g., Shapiro and Campillo, 2004) and to use them for the crustal seismic tomography (e.g., Shapiro et al., 2005; Sabra et al. 2005; Yang et al., 2007). Extracting waves propagating between pairs of stations from noise correlations also provides us with a possibility to have repeatable measurements of travel times and to detect their changes related to the perturbation of the media. This gives rise to noise based seismic monitoring (e.g., Sens-Schönfelder and Wegler, 2006; Brenguier et al., 2008a,b).

The noise based Green function reconstruction implies, however, some strong hypothesis about the noise modal composition. A perfect reconstruction can be achieved for an ideally random and equipartitioned wavefield (Sánchez-Sesma and Campillo, 2006). In a case of seismic noise it would imply that it sources should be distributed randomly and homogeneously in volume. This is obviously not the case for the real ambient seismic noise within the Earth. First, most of its sources are located on the surface resulting in stronger presence of fundamental mode surface waves and their relatively easy reconstruction from inter-station cross-correlations. Second, the distribution of noise sources is not perfectly random and homogeneous. Background seismic oscillations are mostly generated by the forcing from oceanic gravity and infragravity waves. The interaction between these oceanic wavescand the solid Earth is governed by a complex non-linear mechanism (e.g., Longuet-Higgins, 1950) and, as a result, the noise excitation depends on many factors such as the intensity of the oceanic waves but also the intensity of their interference as well the seafloor topography (e.g., Kedar et al., 2007). Overall, the generation of seismic noise is strongly modulated by strong oceanic storms and therefore, has a clear seasonal and non-random pattern. Distribution of noise sources homogenizes when considered over long times (more than one year). The homogenization and randomization of the noise wavefield is also enhanced by the scattering of the seismic waves on the small-scale heterogeneity within the Earth. Also, because of the stationary phase principle, a cross-correlations of the noise recorded by two receivers is dominated by contribution from sources located in vicinity of the line connecting these receivers. Therefore, even without a perfectly homogeneous distribution, a presence of sufficient amount of favorably located noise sources results in relatively high quality reconstruction of fundamental mode surface waves. As a consequence, reconstructing surface waves from correlations of seismic noise and measuring their dispersion curves works rather well. However, further improving the accuracy of the noise based measurements needed to develop new high-resolution imaging and monitoring methods requires better understanding of the noise modal content and its evolution in space and time. Taking into account realistic distribution of noise sources is also necessary to be able to apply waveform inversion approaches to the noise correlations.

References
Brenguier, F., N. Shapiro, M. Campillo, V. Ferrazzini, Z. Duputel, O. Coutant, and A. Nercessian (2008a). Towards forcasting volcanic eruptions using seismic noise. Nature Geosciences. doi :10.1.1038/ngeo104.
Brenguier, F., M. Campillo, C. Hadziioannou, N. M. Shapiro, R. M. Nadeau, and L. E. (2008b). Postseismic Relaxation Along the San Andreas Fault at Parkﬁeld from Continuous Seismological Observations. Science 321, 1478.
Derode, A., E. Larose, M. Tanter, J. de Rosny, A. Tourin, M. Campillo, and M. Fink (2003). Recovering the Green’s function from ﬁeld-ﬁeld correlations in an open scattering medium (L). Acoustical Society of America Journal 113, 2973–2976.
Gouédard, P., L. Stehly, F. Brenguier, M. Campillo, Y. Colin de Verdière, E. Larose, L. Margerin, P. Roux, F. J. Sánchez-Sesma, N. M. Shapiro, and R. L. Weaver (2008). Cross-correlation of random ﬁelds: mathematical approach and applications. Geophysical Prospecting 56.
Kedar, S., M. Longuet-Higgins, F. Webb, N. Graham, R. Clayton, and C. Jones (2007). The origin of deep ocean microseisms in the North Atlantic Ocean. Proc. R. Soc. A. doi:10.1.1098/rspa2007.0277. Longuet-Higgins, M. (1950). A theory of the origin of microseisms. Phil. Trans. R. Soc. A 243, 1 – 35.
Sabra, K., P. Gerstoft, P. Roux, W. Kuperman, and M. Fehler (2005, April). Surface wave tomography from seismic ambient noise in Southern California. Acoustical Society of America Journal 117, 2431–2432.
Sánchez-Sesma F J and Campillo M. (2006). Retrieval of the Green function from cross-correlation: The canonical elastic problem. Bull. Seism. Soc. Am. 96: 1182-1191. Sens-Schönfelder, C. and U. Wegler (2006). Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia. Geophysical Research Letters 33, 21302.
Shapiro, N. M. and M. Campillo (2004). Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise. Geophysical Research Letters 31, 7614.
Shapiro, N., M. Campillo, L. Stehly, and M. Ritzwoller (2005). High-resolution surface wave tomography from ambient seismic noise. Science 307, 1615–1618.
Toksöz, M. (1968). Microseisms : Mode structure and sources. Science 159, 872–873.
Yang, Y., M. Ritzwoller, A. Levshin, and N. Shapiro (2007). Ambiant noise Rayleigh tomography across Europe. Geophys. J. Int. 168, 259 –274.

Top

Seismic tomography of Eurasia
Tsvetcova Т., Shumlyanska L., Zaiets L., Bugaienko I. Institute of geophysics of Academy of science of Ukraine, av.Palladina 32, Kyiv tsvet@igph.kiev.ua

The methods of seismic tomography has two tendencies. The first one deals with definition of the relative corrections to the given 1-D reference velocity model (Romanov V., Lavrentev М. [1], Aki K., Christoffersson A., Husebye E. [2]). The second one is based on the idea of the common field of times of wave’s arrival in which the velocity is independent from 1-D reference model are defined. This tendency is represented by methods of Taylor's approximation developed by Geyko V., continuation of the seismic fields by Pillipenko V. [3], method of Goldin S. [4].

3D mantle P-velocity model of Eurasia is obtained by a method of Taylor's approximation developed by Geyko V. [5, 6], the numerical solution is described in papers [7, 8]. The general approach to the solution of a main task of seismic tomography by this method is following:
1) Construction of the common field of an average point of times of arrival of waves at station of supervision;
2) Construction of travel-time curve (TTC) of the refracted wave in the format of an average point;
3) Reversals of the TTC - sections.
The method of Taylor's approximation of the solution of
a seismic tomography task has a number of advantages compared to the classical linearization:
1) It gives essential advantage in accuracy of approximation of nonlinearity, i.e. the received 1-D velocity dependence for each section of the common field of an average point has the same meaning, as 1-D reference model, which is used for all area in a method of classical linearization;
2) The approximation is still correct at weaker restrictions to the function of velocity;
3) It enables resolution of a task of waveguides;
4) It is possible to carry out the common manipulation of TTC of the refracted and reflected wave (from subhorizontal borders);
5) It does not require a choice of 1-D reference of model;
6) The method of Taylor's approximation and its numerical realization is corrected by Tihonov;
7) It considerably reduces dimension of a task of the numerical manipulation;
8) It is equally good for representation of the decision, both in the Cartesian system, and in spherical.

For the construction of 3D mantle P-velocity model the times of arrival of Р-waves published in the ISC bulletins for the period from 1964 to 2004 were used. We used earthquakes with following restrictions: depth of the hypocenter should be less than 50 km; magnitude mb  4.5; number of stations, recorded earthquake 300. The database answers 20000000 seismic lines. The solution is submitted on a grid 0.50.525km. An error of definition of velocity for TTC-sections is 0.015 km/s.

Fig.1. shows the received sections of the common field in a format of an average point of the times of arrival. The given circuit demonstrates incompleteness and non-uniformity of distribution of the data. Are in the best way covered region of Mediterranean, Western Europe, Alps-Himalaya mobile belt (sizes of the area of sections from 0.250.25), the regions Central and East Europe (sizes of areas of sections up to 3.03.0) are worse. Worse all the regions of Siberia and China are covered.

Fig. 1. Sections of the generalized field of times of arrival of Р-waves for: a - Europe; b - Asia.

More than 1100 TTC were used in a form of the sections of the common field in format of an average point that allowed us to reveal the basic regularities of a velocity structure of a mantle of Eurasia. Incompleteness of the data resulted in non-uniform investigation of a mantle on depth, in a southern half of researched area the model is constructed up to depth of 2500 km while in north it was developed only to 850 km.

The received results are presented in horizontal, longitudal and latitudinal sections. The horizontal ones are presented in true speeds through 25 km, vertical - through 1 on all area of research in deviations concerning average velocity (δVp=Vav-Vp), received for mantle of Eurasia. The vertical sections are submitted on depth from 850 km to the north 70 n.lat., in an interval 60 – 70 n.lat. - up to 1100 km, in an interval 50 - 60 n.lat. - 1700 km, to the south 50 n.lat. - up to 2500 km. All 455 sections are received.

Fig. 2. shows correspondence between theoretical and natural TTC.

Fig. 3. Correspondence between theoretical and experimental TTC: a - one from a TTC from the Black Sea, another one represents the TTC of the Mediterranean area.

The analysis of the obtained travel-time curves for the territory under Eurasia shows the following points of break: 18° -25° – correspond to the upper part of the mantle and transition zone; 62° -65° – corresponds to the border of the division between zone ІІ and a bottom part of the middle mantle; 84° - 85° - allocates subdivisions in the lower mantle. The obtained breaks of the TTC correspond to the accepted borders in a mantle [9}. The borders change depending on region.

The obtained velocity model of the mantle of Eurasia has the following features:
- the upper mantle and partially zone of division І is essentially non-uniform, average and bottom mantle quasihomogeneous with allocation of more non-uniform part on depth of 1600-1700 km (the zone of division - ІІ);
- Two types of Precambrian platforms with the various mantle velocity characteristics were recognized. The first type includes East-European, northern part of the Siberian, the Indian platforms, which are characterized by the greatest thickness of a high-velocity layer of seismological lithosphere. Average thickness of the seismological lithosphere of the Eastern Europe and northern part of the Siberian platform is 350 km, and 200 km – for the Indian platform. Another characteristic is a low-velocity nature of a transition zone. Second type of platforms includes South-Chinese and North-Chinese platforms, which have moderate thickness (up to 100 km) of a high-velocity layer of seismological lithosphere and low-velocity characteristic of the mantle down to 450 km; transition zone under the this type of platforms represented by a high-velocity layer.
- mantle of Phanerozoic structures have inverted nature in comparison to Precambrian platforms, i.e. upper mantle is of low-velocity nature while transition zone is high-velocity. By these features the Eurasians belt was allocated;
- the inclined position of a high-velocity layer of seismological lithosphere observed that define basically borders of the Eurasians belt (fig. 3.)
- on the background of a quasihomogeneous middle mantle and zone of II division velocity anomalies correlated to the cold Deccan plume can be seen. These anomalies deviate from the average high-velocity background by δVp0.05 km/s. The anomalies on an average high-velocity background of the middle and lower mantle correspond to known plumes - Afar, Siam, etc. The sharp low-velocity anomaly is allocated under the Black sea in the lower mantle. Down to the depth of 1500 km inclined high-velocity anomalies corresponding to "slabs" are allocated.

Fig.3. Borders of the Eurasians belt.

1. Lavrentev М., Romanov I. About three linearized return tasks for the hyperbolic equations.// Report of academy of science of USSR, 1966,171, P.1279-1281.
2. Aki. V., Christoffersson,A., Husebye, E. Determination of the three-dimensional seismic structure of the lithosphere // J.Geoph.Res, 1977,82, P.277-296.
3. Pillipenko V. Numerical method of fields of times for construction of seismic borders / Inverse kinematical tasks of explosion seismology, M:” Nauka”, 1979, 232 p.
4. Goldin S. Inverse problems ray seismic tomography// Geology and geophysics, 1997, 38, №5, P.981-998.
5. Geyko V. Taylor′s approximation of the wave equation and equations eikonal in inverse seismic problems // Geophys. journal, 1997, 19, №3, P.48-68.
6. Geyko V. A general theory of the seismic travel-time tomography // Geophys. journal, 2004, 26, №1, P.3-32.
7. Geyko V., Tsvetcova T. About uniqueness of the solution one-dimensional inverse kinematical seismic problem // Geophys. journal, 1989, 11, №6, P.61-67.
8. Tsvetcova Т. Continuation of the common field of times downwards and inverse kinematical seismic problem // Report of Academy of science of Ukraine, 1996, №1, P.79-85.
9. Puscharovskyy J. Seismic tomography, tectonics and deep geodynamics // Report of Academy of science of Russia, 1998, 360, №4, P.518-522.

Top

Processes of percolation in rocks
Sidorova I.P. Institute of Geology and Geophysics of Academy of Sciences of Uzbekistan, Tashkent, irina_sid@mail.ru

At modeling processes of transportation of heat from Earth’s depths the dominating attention of researchers occupies the convection mechanism of heat transfer. At the same time the conduction plays not less important role and is more complex phenomenon for understanding, because the rocks, through which there is transferred the heat, are complex multicomponent natural ensembles.
In this connection in laboratory of Regional and Applied geophysics of Institute of Geology and Geophysics of the Academy of Sciences of Uzbekistan the computer program «PercolationCheck» (the patent №DGU 01170), intended for calculation of percolation (conductivity) process in rocks on the basis of raster pictures of distribution of elements in the sample, received on microanalyzer JXA-8800R JEOL has been developed. The scope of the program - research of materials - chemistry, physics of natural combinations. The program determines presence cells of percolation configurations and visually displays on the screen of computer the available in the investigated sample the clusters, their capacity, presence of one or more connecting clusters, testifying about presence's percolation (conductivity).
To development preceded long-term (more than 20 years) laboratory studying of thermophysical properties of rocks. For this period data about heat conductivity, thermal capacities from superdeep Muruntau borehole SG-10 [1, 2, 3], and also by ore deposits and the areas of the Western Uzbekistan adjoining to them are in details studied and generalized. Scientific researches have been directed on the analysis of the reasons of occurrence thermophysical anomalies in the certain rocks on the basis of physical and mathematical modelling processes of conductivity in heterogeneous environments – in rocks. The research problem of thermophysical anomalies in rocks were solved in four stages:
1) thermophysical experiments on hardware complex ITEM-1М, IТ-λ-400 and IТ-с-400; 2) the mineralogical analysis of samples on microanalyzer JXA-8800R JEOL; 3) development of physical and mathematical model of heat conductivity of rocks; 4) revealing of anomalous geoobjects on the basis of application of algorithm of the program «Percolationcheck» on experimental data.
As a result of investigations of the first stage have been revealed anomalous (on heat-conducting properties) groups of rocks from collection by SG-10. Samples investigated by us macroscopic were identical, however at detailed studying at the second stage was elucidated, that in different samples distribution of components on volume absolutely variously. Researches on the allocated ensembles at III stage were spent on polished sections on microanalyzer JEOL JXA - 8800R, time of the analysis of the chemical compound of separate mineral phases on EDS Link ISIS-300 (Oxford, the Great Britain) made from 1 till 5 minutes. As a result of the lead experiments on IY stage a generalizing stage basic elements which distribution in the sample essentially influences conductivity of heat are allocated. From the physical point of view they form in samples infinite clusters is fixed by the program «Percolationcheck», basic which working modules are shown on figure 1: 1) module «Geodata» provides input of the raster information; 2) in module «Calibration» there is an allocation of ranges; 3) in module «Сluster» are allocated conducting clusters; 4) module «Stat» conduct statistical processing; 5) module «Result» classifies samples on conductivity properties.

Figure 1. The scheme of realization of the program «Percolationcheck»
     a) raster picture of distribution of element carbon, received on the microanalyzer;
     b) visual display infinite cluster;
     c) the basic working modules of the program «Percolationcheck»

At computer modeling process of percolation we used model «cell’s percolation», in which each cell can be in two conditions – «occupied» or «empty». Thus each cell is engaged with probability р irrespective of condition of the next cells. The occupied cells either are isolated from each other, or form the groups consisting of the nearest neighbors. Cluster is defined as group of the occupied cells of the lattice connected with nearest neighbors on the side of cell. The program «Percolationcheck» generates cell’s percolation configuration and displays its on monitor of computer for visualization of representation clusters. After marks of clusters we receive need geometrical characteristics. For reception of the cluster's size is summarized on correct labels, thus, as a result of calculations on the algorithm incorporated in the program, the estimation of prognostic occurrence percolation processes (anomalous conductivity) on raster pictures of distribution of elements in samples of rocks is made.

References
1. Sidorova I.P. Thermophysical complexes in the section of the Muruntau deep borehole SG-10 (Central Kyzylkums, Western Uzbekistan), Geologiya va mineral resurslar, 2000, №1, рр.33-38 [in Russian].
2. Sidorova I.P. Thermophysical anomalies in rocks, uncovered by deep borehole SG-10 (Central Kyzylkums, Western Uzbekistan).- Abstract /Report/-XXIII General Assembly of International Union of Geodesy and Geophysics (IUGG), Sapporo, Japan, 2003.
3. Sidorova I.P. Thermophysical anomalies in rocks, uncovered by deep borehole SG-10 /Central Kyzylkums, Western Uzbekistan/. Abstract /Report/ 32Pnd P International Geolgical Congress Florence, Italy, 2004.
4. Efros A.L. The Physics and Geometry of the Chaos// M.: Science, 1982, 260p.
5. Birch F., Clark H. The Thermal Conductivity of Rocks and its Dependence upon Temperature and Composition// Amer.J.Sci.,1940, V.238, pp.529-558.

Top

Initiation of the modern style of subduction in the Precambrian: insights from numerical experiments
E. Sizova¹, T. Gerya^1,2, B. Kaus¹, M. Brown³, L. Perchuk^4,5
1. Institute of Geophysics, ETH-Zurich, Sonnegstrasse, 5, 8092 Zurich, Switzerland (sizova@erdw.ethz.ch),
2. Geology Department, Moscow State University, 119199 Moscow, Russia,
3. Department of Geology, University of Maryland, College Park, MD, 20742, USA,
4. Department of Petrology, Moscow State University, 119199 Moscow, Russia,
5. Visiting Professor of Faculty of Science, University of Johannesburg, South Africa

Plate tectonics is a self-organizing global system driven by the negative buoyancy of the thermal boundary layer resulting in subduction. Although the signature of plate tectonics is recognized with some confidence in the Phanerozoic geological record on continents, evidence for plate tectonics is less certain further back in time. To improve our understanding of plate tectonics on the Earth during the Precambrian, we have to combine knowledge derived from the geological record with results from realistic numerical modeling.

In a series of experiments using a 2D petrological–thermomechanical numerical model of oceanic-continental subduction we have systematically investigated the dependence of tectono-metamorphic and magmatic regimes at an active continental margin on upper-mantle temperature, crustal radiogenic heat production, degree of lithospheric weakening as well as other physical parameters. The model includes spontaneous slab bending, dehydration of subducted crust, aqueous fluid transport, mantle wedge melting, and melt extraction from the mantle resulting in crustal growth. We have identified a first-order transition from a “no-subduction” tectonic regime through a “pre-subduction” tectonic regime to the modern style of subduction. The first transition is gradual and occurs at upper-mantle temperatures between 250-200 K above present-day values, whereas the second transition is more abrupt and occurs at 175-160 K. The link between geological observations and model results suggests that the transition to the modern plate tectonics regime might have occurred during the Mesoarchean–Neoarchean time (ca. 3.2–2.5 Ga). In the case of a “pre-subduction” tectonic regime (upper-mantle temperature 175–250 K above present) the plates are weakened by intense percolation of melts derived from the underlying hot melt-bearing sub-lithospheric mantle. In such cases, convergence does not produce self-sustaining one-sided subduction, but rather results in shallow underthrusting of the oceanic plate under the continental plate. A further increase in the upper-mantle temperature (> 250 K above present) induces a transition to a “no-subduction” regime where horizontal movements of small deformable plate fragments are accommodated by internal strain and even under imposed convergence shallow underthrusts do not form.

To better understand the underlying physics of these models we performed an additional series of experiments using similar 2D petrological–thermomechanical numerical model but without hydration, melting and extraction procedures. In these models, we have obtained a similarly abrupt transition from the modern style of subduction to the “no-subduction” regime at the upper-mantle temperature 160-180 K above the present-day values. This temperature is approximately the same as determined in the first set of experiments. The “no-subduction” regime is characterized by ‘dripping-off’ of the plate tips, most likely because of the small effective viscosity constrast between subducting slab and surrounding mantle. Indeed we do not observe a transitional “pre-subduction” tectonic regime with underthrusting of the oceanic plate in these sets of models. This implies critical role of rheological weakening by sublithospheric melts in defining how transition between ancient “no-subduction” stage and modern plate tectonic regime occurred in the Earth’s history.

Top

Three dimensional numerical modeling of lithospheric dynamics coupled with mantle convection
Stephan V. Sobolev, Anton Popov and Bernhard Steinberger GFZ German Research Centre for Geosciences (GFZ), Potsdam (stephan@gfz-potsdam.de)

The convection in deep Earth is connected to the surface through the heterogeneous and rheologically complex lithosphere and asthenosphere. Their rheology is usually oversimplified in global geodynamic models, which limits possibility to model plate tectonics. In this work, we use our experience in modeling of lithospheric deformations with complex rheology to develop more realistic models of coupled global mantle convection and 3D lithospheric dynamics.

We employ the newly developed 3D thermomechanical finite element numerical technique to model a 300 km thick upper layer of the Earth, coupled with the convecting mantle. The present day temperature distribution and crustal structure within the layer are taken from existing models. We also assume that the upper layer is composed from rocks with non-linear temperature- and stress-dependent visco-elastic rheology, corresponding to the dry or wet olivine (mantle) or naturally wet plagioclase (crust), combined with Mohr-Coulomb frictional plasticity. Plate boundaries are represented by the narrow zones of elasto-visco-plastic rheology with much lower frictional strength than within the plates. The mantle below the 300 km depth is modeled using Hager and O’Connell’s mantle flow spectral modeling technique with present day density and viscosity distribution based on either interpretation of global seismic tomography or history of subduction. The upper layer and mantle modeling domains are coupled by iteratively achieved continuity of tractions and velocities at 300 km depth.

Here we will show modeling results focusing on the effect on the plate velocities of lateral viscosity variations in the mantle and of the frictional strength at plate boundaries. Particularly, we consider the effect of the viscosity of the asthenosphere on plate velocities and net rotation of the lithosphere. Modeling shows that deep convection generates plate tectonic-like velocity pattern only when effective friction at plate boundaries becomes less than 0.05-0.1. Both magnitudes and directions of plate velocities are reproduced very well at friction in subduction zones of less than 0.03-0.05 and friction at other plate boundaries of 0.05-0.1. The best fit of the observed velocities and net rotation of the lithosphere is obtained assuming that asthenosphere is significantly more “wet” than lithosphere. We also show that lateral viscosity variations in the mantle deeper than 300 km do not affect plate velocities, while such variations in the upper 300 km have very significant influence on magnitudes of plate velocities.

Top

Plume formation in strongly temperature-dependent viscosity fluids and the origin of the Martian hemispheric dichotomy
Slava Solomatov and Yun Ke Department of Earth and Planetary Sciences, Washington University, St. Louis, Missouri 63105, U.S.A

One of the most prominent features of Mars is the hemispheric dichotomy. The Martian surface consists of a heavily cratered elevated southern hemisphere and a resurfaced depressed northern hemisphere. The dichotomy seems to have formed very early in the history of the planet (Frey, 2006). Several hypotheses have been proposed to explain the origin of the dichotomy including exogenic hypotheses (Wilhelms and Squyres, 1984; Frey and Schultz, 1988; Reese and Solomatov, 2007; Andrews-Hanna et al. 2008; Marinova et al., 2008; Nimmo et al., 2008) and endogenic hypotheses (Wise, 1979; McGill and Dimitriou, 1990; Sleep, 1994; Zhong and Zuber. 2001; Elkins-Tanton et al., 2003, 2005; Roberts and Zhong 2006). However, the issue remains poorly understood. One of the main difficulties is the complexity of physical processes on early Mars. We focus on one particular process, plume formation and propose that the dichotomy can be caused by a transient superplume produced by a hot Martian core (Ke and Solomatov, 2004, 2006, 2009).

The number of plumes on early Mars would be quite large if the Martian mantle were a constant viscosity fluid. However, the viscosity of rocks varies by many orders of magnitude. Fluid dynamics of plume formation in such fluids depends on various parameters of the system. We investigate plume formation in strongly temperature-dependent viscosity fluids using theoretical analysis, parameterized convection calculations and finite element calculations in fully three-dimensional, spherical shell geometry by CitcomS (Zhong et al., 2000). We find that when the viscosity contrast across the thermal boundary layer at the base of the mantle is sufficiently high, a large single plume (.superplume.) forms. Plume formation involves several processes including crystallization and growth of the thermal boundary layer (TBL) at the base of the mantle, small-scale convection within the TBL and large-scale instability of the TBL. The superplume can easily satisfy the timing constraints on the formation of the dichotomy. Also, the core cooling is sufficiently rapid to induce convection inside the core and allow the operation of the magnetic dynamo. Depending on the model parameters, the magnetic field exists for millions to hundreds of millions of years after planetary formation, consistent with observations.

An important condition for superplume formation is that the early Martian core was hotter than its mantle at least by several hundred degrees Kelvin. This is supported by the following arguments. (i) On the present-day Earth, the temperature contrast between the lower mantle and the core is estimated to be 1000 to 2000 K (Boehler, 1996, 2000; Williams, 1998) and was probably even higher in the past as the core temperature must have decreased more than the mantle temperature during planetary evolution (Stevenson et al., 1983). (ii) Numerical simulations of giant impacts show that during collisions of large differentiated bodies with the proto-Earth the material that was heated most was the iron core of the impactor (Canup, 2004). After the iron from the impactor's core segregated into the Earth's core, the temperature contrast between the core and the mantle reached several thousand degrees Kelvin. For a smaller planet like Mars the core-mantle temperature difference after a giant impact may be high as well, of the order of 1000 K (Ke and Solomatov, 2006). (iii) A substantial fraction of the gravitational potential energy of core formation can, in principle, be retained within the segregated iron (Stevenson, 2001). Although there are many factors controlling the energy partitioning between the mantle and the core, how liquid iron segregated might be the most important one. If iron moves in the form of blobs, then the gravitational potential energy is deposited largely withing the mantle. Very little temperature difference is expected in this case. However, numerical simulations show that magma migration (Wiggins and Spiegelman, 1995; Connolly and Podladchikov, 2007) as well as segregation of liquid iron (Golabek et al., 2008) results in the formation of melt channels. If liquid iron does segregate through channels then the viscous dissipation occurs largely within the channels and gravitational potential energy dominantly goes into the core (Ke and Solomatov, 2009). (iv) An overturn of the Martian mantle would adiabatically cool the lower part of the mantle and increase the core-mantle temperature contrast even further. According to Elkins-Tanton et al. (2005), the temperature contrast between the lower mantle and the core due to mantle overturn can reach nearly 1000 K. A combination of the effects described above may generate the temperature contrast between the core and the lower part of the Martian mantle from several hundred to several thousand degrees Kelvin. According to our theory, this is sufficiently high to generate a superplume in the early history of Mars and explain the formation of the dichotomy.

Top

Fluid flow dynamic in transient growing regime of layered visco-elastic saturated sediments.
Suetnova E.I. Institute of the Physics of the Earth, RAS. Moscow, Russia. elena _suetnova@mail.ru

Knowledge of the dynamic of fluid flow controlling the distribution and the transport of liquid in the Earth crust is necessary to quantitatively model many geological processes. An amount of liquid in rock modifies its physical properties, like its electrical and thermal properties and seismic wave propagation. Fluid flow plays a key role in formation of sedimentary basin feel during geological time of sedimentary basin formation. The compaction of accumulated sediments is main process, which lead to flow of fluid through sediment during basin fill forming. Sediment compaction and associated fluid flow in sediments have been studied under a various approaches since the middle of 20 century. The recent results in this area include the study of the viscous sediment compaction [1-4], the viscous compaction wave development in sediments [5, 6], visco-elastic sediment compaction [7-8]. In the approach [8], it has been shown the dependence of fluid flow dynamic upon sedimentation rate or velocity of sediment peel growing during basin formation. Sediments have a complex rheology, in which the balance between the elastic and viscous contribution to the deformation is time-dependent [8]. The fluid flow in saturated sediments depends on effective pressure gradient and sediment permeability. Because permeability nonlinearly depends upon porosity and effective pressure gradient are interrelated with deformation of sediment matrix under the loading, the dynamic of fluid is controlled by the behavior of deformable matrix. Mathematical formulation of physical problem of low-viscosity fluid flow in a more or less deformable visco-elastic matrix with moving boundary consist of system of nonlinear partial differential equations with representative boundary conditions [8]. Analysis of dimensions of variables and parameters (dimensional analysis) allow formulating the similarity numbers, which determined the regime of porosity and fluid flow evolution during the time of sediment growing for given rheological and physical property of fluid and sediments and sedimentation rate. The fluids flow is stable in the case of uniform sedimentation and impermeable lower boundary. Unfortunately, uniform sedimentation process unlikely corresponds to a real process of sediment accumulation. Geophysical data show that the sediments type varies through sediments column. To study flow dynamic and visco-elastic compaction in the case when different type of sediments have been accumulated during different time intervals of basin history, the system of governing equations have been modified [9] by introducing the time intervals of given types of sedimentations during sedimentation history. Model calculations show that the pressure seals and pressure cell are formed in dependence of sedimentation rate and permeability or viscosity contrast. The time of development of such peculiarity depend on sedimentation rate and rheological properties of sediments, the amplitude of pressure and fluid flow variations depend on permeability or viscosity contrast and the time of sediment fill growing. The characteristics of the flow anomaly development depend upon the sedimentation rate. Fig.1 illustrates two different cases of development of fluid pressure anomalies for sedimentation rate 10-11m/s (all carve are dimensionless). At the fig. 1a it is shown the evolution of fluid pressure in the case when relatively higher-viscosity sediments accumulate in the basin within the time interval over which the sedimentary thickness increases from 0.15 to 0.2 of its final value. At the fig. 1b it is shown the evolution of the depth distributions of fluid pressure in the case when lower-viscosity sediments accumulate in the basin within the time interval over which the sedimentary thickness increases from 0.35 to 0.37 of its final value. Some spatial scenario of sedimentation history may leads to generations of flow instability and hydro-fracturing (fig.1 b). Fig.2 illustrates the dependence of transient fluid flow in growing layered saturated sediments upon sedimentation rate. At the fig. 2a it is shown the evolution of the depth distributions of fluid pressure in the case when relatively lower permeability sediments accumulate in the basin within the time interval over which the sedimentary thickness increases from 0.17 to 0.5 of its final value. At the fig. 2b it is shown the evolution of the depth distributions of fluid pressure in the same case when lower- permeability sediments accumulate in the basin within the time interval over which the sedimentary thickness increases from 0.17 to 0.5 of its final value and sedimentation rate increases up to 10-10m/s (actual final sediment thickness is equal in both cases). Thus, the similarity numbers determines the dynamic of fluid flow disturbances during accumulation of layered sediments.

References.
1. McKenzie, D. The compaction of igneous and sedimentary rocks // J. Geol. Soc. London. 1987. V. 144. P. 299–307.
2. Birchwood, R.A. and Turcotte, D.L. A unified approach to geopressuring, low-permeability zone formation, and secondary porosity generation in sedimentary basins // J. Geophys.Res.1994. V. 99. P. 20 051–20 058.
3. Fowler A.C., Yang X. Pressure solution and viscous compaction in sedimentary basin // J. Geophys. Res. 1999. V. 104. P. 12989-12997.
4. Suetnova E. I. and Chernyavskii V. M. Free fluid Flow through a visco-deformable porous medium with a moving boundary // Fluid Dynamics. 2001. V. 36. N. 1. P. 121–129.
5. Karakin A. V. and Suetnova E. I. Simulation of a flow wave in a porous saturated visco-deformable medium // Izv. Akad. Nauk SSSR, Fizika Zemli. 1989. N. 3. P. 94–99.
6. Connolly J.A.D., Podladchikov Yu.Yu. Temperature-dependent viscoelastic compaction and compartmentalization in sedimentary basins // Tectonophysics. 2000. 324. P. 137-168.
7. Yang, X.-S. Nonlinear viscoelastic compaction in sedimentary basins // Nonlinear Proc. Geophys. 2000. 7. P. 1–7.
8. Suetnova E.I. and Vasseur G. 1-D Modelling rock compaction in sedimentary basin using viscoelastic rheology // Earth Planet. Sci. Lett. 2000. V. 178. P. 373–383.
9. Suetnova E.I. Viscoelastic compaction of heterogeneous sediments // Izvestiya, Physics of the Solid Earth. V. 39. N. 1. 2003. P. 71–77.

Top

3-D spherical models of coupled mantle thermo-chemical evolution, plate tectonics, magmatism and core evolution incorporating self-consistently calculated mineralogy
Paul Tackley¹, Takashi Nakagawa¹, Frédéric Deschamps¹, James Connolly³
1. Institute für Geophysik, ETH Zürich, Sonneggstrasse 5, CH-8092, Zürich (ptackley@ethz.ch)
3. Insitute für Mineralogie und Petrographie, ETH Zürich, Clausiusstrasse 25, CH-8092, Zürich

High pressure and temperature experiments and calculations of the properties of mantle minerals show that many different mineral phases exist as a function of pressure, temperature and composition [e.g.Irifune and Ringwood, 1987], and that these have a first-order influence on properties such as density, which has a large effect on the dynamics, and elastic moduli, which influence seismic velocity. Numerical models of thermo-chemical mantle convection have typically used a simple approximation to treat these complex variations in material properties, such as the extended Boussinesq approximation. Some numerical models have attempted to implement multiple, compositiondependent phases into thermo-chemical mantle convection [e.g., Tackley and Xie, 2003] and to calculate seismic anomalies from mantle convection simulations based on polynominal fitting for temperature, composition and mineral phase [Nakagawa and Tackley, 2006]. However, their linearised treatments are still approximations and may not adequately represent properties including effect of composition on phase transitions. In order to get closer to a realistic mineralogy, we here calculate composition-dependent mineral assemblages and their physical properties using the code PERPLEX, which minimizes free energy for a given combination of oxides as a function of temperature and pressure [Connolly, 2005], and use this in a numerical model of thermo-chemical mantle convection in a three-dimensional spherical shell, to calculate three-dimensionally-varying physical proporties. In this presentation we compare the results obtained with this new, self-consistentlycalculated treatment, with results using the old, approximate treatment, focusing particularly on thermo-chemical-phase structures and seismic anomalies in the CMB region and the transition zone [Nakagawa and Tackley, 2009]. The numerical models treat the evolution of a planet over billions of years, including self-consistent plate tectonics arising from plastic yielding, melting-induced differentiation, and a parameterised model of core evolution based on heat extracted by mantle convection. This approach is also being used to study Mars, Venus, Mercury and super-Earths.

Figure 1. Simulation with PERPLEX-calculated mineralogy and physical properties, from Nakagawa & Tackley, 2009.

References
Connolly, J. A. D., (2005), Earth Planet. Sci. Lett., 236, 524-541.
Irifune, T., and Ringwood, A. E., (1987), Earth. Planet. Sci. Lett., 86, 365-376.
Nakagawa T., and Tackley, P.J., (2006), Geophys. Res. Lett., 33, L12S11, doi:10.1029/2006GL025719.
Nakagawa, T., and Tackley, P. J., (2009), Geochem. Geophys. Geosys., 10, Q03004, doi:10.1029/2008GC002280.
Tackley, P. J., and Xie, S., (2003), Proc. 2nd MIT conference.

Top

Coupling of ¨Mechanical¨ and ¨Chemical¨ Compaction Using Thermodynamic Tools
Evgeniy V. Tantserev¹, Yuri Y. Podladchikov², Christophe Y. Galerne³
1. Faculty of Engineering Science and Technology, NTNU, Norway
2. Physics of Geological Processes, University of Oslo, Norway
3. Steinmann Instüt, Universität Bonn, Germany

Mechanical compaction is usually treated under isothermal (Wilmanski 2006) or isoentropic (Mckenzie 1984) simplifying assumptions. The case of joint consideration of both mechanical compaction and reactive porosity alteration requires somewhat greater than usual care about thermodynamic consistency. There is no “conservation of porosity law”; there are no equilibrium conditions for porosity similar to absence of thermal gradients or continuity of stresses and chemical potentials. Porosity fields may spontaneously develop jumps that do not have to disappear at equilibrium (Hills and Roberts 1988) and may generate ever growing waves out of minuscule perturbations (Connolly and Podladchikov 2007).

Modeling of multi-component multi-phase porous systems is fundamental to the study of geological processes. Recipes of continuum mechanics and thermodynamics are employed here in order to derive a closed system of equations for coupling of different processes. We start from a set of balance equations for mass, momentum, energy, and entropy written for a general multi-phase multi-component system. Based on a local thermodynamic equilibrium assumption and non-negativity of entropy production, we choose thermodynamically admissible expressions for fluxes and sources in the set of balance equations. To simplify treatment of the chemical reactions, we assume local chemical equilibrium reached in each point due to the high temperature of the system and slow rate of percolation as it is assumed in an ideal chromatography theory, e.g. (Helfferich and Whitley 1996). Our goal is to derive a thermodynamically admissible closed system of equations (CSE) describing the coupling of “mechanical” and “chemical” compaction, e.g. (Cass T. Miller and Gray 2008, Helfferich and Whitley 1996, Spiegelman et al. 2001). We are restricting ourselves to a minimum set of the most essential processes. The emphasis is not on generality or universality of the final system of equations but on its transparency and thermodynamic consistency. For the mechanical interactions, we include pore compressibility and viscosity to account for high temperature stress relaxation and Darcian flow of the porous fluid. Due to their complexity, we exclude the capillary effects. For the energy balance, we introduce a specific (per unit mass) internal energy for each phase and assume local thermal equilibrium that requires equal temperatures for all phases. Our CSE includes evolution equation for porosity. We calibrate our poroelastic parameters by exact Gassmann’s relationships, e.g. Gassmann (1951), and obtain closing poroelastic rheological relationships.

Our multi-component, compressible, and non-isothermal CSE unify most of the effects considered in formulations of Šrámek al., (2007), Spiegelman et al., (2001) and Aharonov et al., (1995) with the exception of surface tension and non-equilibrium chemical reactions.

Aharonov, E., Whitehead, J.A., Kelemen, P.B. & Spiegelman, M. 1995. Channeling Instability of Upwelling Melt in the Mantle. Journal of Geophysical Research-Solid Earth 100, 20433-20450.
Cass T. Miller & Gray, W.G. 2008. Thermodynamically constrained averaging theory approach for modeling flow and transport phenomena in porous medium systems:1. Species transport fundamentals. Advances in Water Resources 31, 577-597.
Connolly, J.A.D. & Podladchikov, Y.Y. 2007. Decompaction weakening and channeling instability in ductile porous media: Implications for asthenospheric melt segregation. Journal of Geophysical Research-Solid Earth 112, -.
Gassmann, F. 1951. Überdie Elastizität poroser Medien. Veirteljahrsschrift der Naturforschenden Gasellschaft in Zürich 96, 1-23.
Helfferich, F.G. & Whitley, R.D. 1996. Non-linear waves in chromatography .2. Wave interference and coherence in multicomponent. Journal of Chromatography A 734, 7-47.
Hills, R.N. & Roberts, P.H. 1988. A Generalized Scheil-Pfann Equation for a Dynamical Theory of a Mushy Zone. International Journal of Non-Linear Mechanics 23, 327-339.
Mckenzie, D. 1984. The Generation and Compaction of Partially Molten Rock. Journal of Petrology 25, 713-765.
Spiegelman, M., Kelemen, P.B. & Aharonov, E. 2001. Causes and consequences of flow organization during melt transport: The reaction infiltration instability in compactible media. Journal of Geophysical Research-Solid Earth 106, 2061-2077.
Šrámek , O., Ricard, Y. & Bercovici, D. 2007. Simultaneous melting and compaction in deformable two-phase media. Geophysical Journal International 168, 964-982.
Wilmanski, K. 2006. A few remarks on Biot's model and linear acoustics of poroelastic saturated materials. Soil Dynamics and Earthquake Engineering 26, 509-536.

Top

EVOLUTION OF MANTLE CONVECTION WITH SELFGENERATING LITHOSPHERIC PLATES AND MOVING CONTINENTS
Trubitsyn V.P., IPE RAS, Moscow

1. Mantle is modeled by two-components viscous liquid with temperature, pressure and stress depended viscosity in a long box 10:1 with grid 800:200 with Ra=107 in expanded Boussinesq approximation by Citcom code (Moresi, Zhong). Floating thick continent is modeled with active markers as buoyant extra high viscous area. For diffusion and dislocation creep mantle convection takes place under motionless thin stagnant lithospheric lid. After Peierls creep was added the lithospheric lid breaks into many quasi-rigid plates. A buoyancy continent (100 km thick and 6000 km long) was put on mantle. Numerical models of mantle convection thermal and mechanical coupled with continent and selfgenerating lithospheric plates show long evolution of this system. Models reveal variations in time of the sizes and number of plates, velocity variations, declination of subducted slabs, back-rolling, change of thermal and stress state of mantle, lithospheric plates and continent, end so on.

2. Models of mantle convection with Newton viscosity with floating continent were calculated using Citcom code and finite differences code(Rykov, Trubitsyn). They simulate main features of Pacific ocean evolution. Model explains extensional and compressional state of mantle between continent and subduction zone, origin and evolution of marginal seas similar to South-East Eurasia boundary, closure of marginal seas similar to S. America west boundary during its motion to the west and N. America encountering with mid-ocean ridge leading to California splitting.

3. 3-D spherical model of mantle convection coupled with 16 rounded continents shows periodical (four times) continental assemblages into supercontinents and then their disperses. 4. 3-D spherical model (58x144x216, Ra=107) simulates continental drift for 150 Ma in future. The initial mantle temperature was calculated from seismic tomography data. Six continents and nine islands with real figures were put on mantle according to their present-day position. Model takes into account thermal and mechanical coupling between continents and mantle and with each others. Calculated heat flow and continent velocities are in accordance with observable data for present-day Earth. Model shows detailed future drift of each continent, future opening Arctic ocean and motions of all continents to assemble them into new supercontinent on south Earth hemisphere.

Top

Discrepancy and seismic isotropy in D''
D.Komatitsch⁽¹⁾, L. Vinnik⁽²⁾ and S. Chevrot⁽³⁾
1. Universite de Pau et des Pays del l´Ádour, France
2. Institute of physic of the Earth, Moscow, Russian Federation, vinnik@ifz.ru
3. Universite Paul Sabatier, Toulouse, France

The D'' layer atop the core-mantle boundary is not only heterogeneous but also anisotropic, at least at some locations. The usually assumed form of anisotropy is hexagonal with a vertical symmetry axis (transverse isotropy). Origins of this anisotropy remained unclear until it was found that at the pressures and temperatures in D'' the silicate perovskite is transformed into the post-perovskite phase with anisotropic crystals that can be preferentially oriented. Anisotropy at the base of the mantle can be an instrument for probing post-perovskite crystal orientations and the related mantle flow. A major seismic technique for measuring anisotropy in D'' is based on observations of Sdiff seismic phases which propagate along the core-mantle boundary. The SVdiff phase with vertical polarization sometimes (but not always) arrives later than the horizontally polarized SHdiff by ~2 seconds, and the difference in the arrival times is interpreted in terms of D'' anisotropy. In this study we examine the strength of evidence for D'' anisotropy and investigate a possibility of observing the SHdiff/SVdiff discrepancy in isotropic D". Using numerical simulations we demonstrate that this discrepancy exists even in a laterally homogeneous D". A delay of SVdiff relative to SHdiff is significant in IASP91 (isotropic Earth reference model) and in other models with a high-velocity D''. The delay is accompanied by a large attenuation of SVdiff with distance. In the models with a low-velocity D'' SVdiff propagates with a low attenuation and with the same velocity as SHdiff or, depending on the vertical S velocity gradient, faster than SHdiff. Propagation of Sdiff in a laterally heterogeneous D'' was simulated with the spectral element technique. In this case SVdiff can be delayed because, unlike SHdiff, SVdiff tends to choose low-S-velocity propagation paths. For the D" models with the S velocity difference of 3% between the hemispheres, there are regions where SHdiff arrives up to ~15 s earlier than SVdiff. The region of early SHdiff arrivals is especially large over the high-velocity hemisphere. If these effects are interpreted in terms of D'' anisotropy, the estimates of anisotropy are biased and may contribute to the popular idea of seismic anisotropy being different in the high- and low-S-velocity regions of D''

Top

Top

Precise stationary gravimetry in studies of the Earth’ inner structure and its dynamics
V.Volkov¹, A.Kopaev²
1- Institute for Physics of the Earth of RAS, Moscow, volkov@ifz.ru
2- Sternberg State Astronomical Institute of Moscow State University, Moscow, gravity3@yandex.ru

In this revue we demonstrate how precise gravity observations on the Earth’ surface give important information about its structure and global geodynamical phenomena.
Relative gravimeters are being used for this purpose starting from 50-s-60-s of the last century, in 70-s special stationary devices have been developed, but the use of all of them is limited by the earth tide investigations in diurnal and semidiurnal bands only. Starting from the 1960-s, IPE RAS participates very actively earth tide studies on the FSU territory as well as abroad. In 1990-s SAI MSU became active in this investigations as well. We will present the main achievements of our institutes in tidal gravimetry.
Currently the most powerful tools for stationary gravity observations are superconducting gravimeter (SG). Their observational capability ranges from seismic frequencies to periods more than 1 year and accuracy varies from 0.1 nGal for quasi-stationary processes to 1 nGal for coherent signals (earth tides and eigenvibrations).
The most intriguing are reports about the observations of eigenvibrations of Earth’ liquid and inner cores, although these results still need confirmation. Study of these phenomena helps to make conclusions about ellipticity of the inner core, viscosity of the liquid core and electromagnetic coupling between liquid core and mantle. In tidal range the most important phenomen is the fine structure of so-called nearly diurnal resonance, that is caused by the interaction between liquid core and mantle. Its investigation helps to make more precise the estimates of ellipticity and quality factor of the liquid core. Comprehensive analysis of long term observations in diurnal and semidiurnal ranges helps to improve the models of the oceanic and loading tides, as well as to estimate the viscosity of the mantle in tidal band and lateral anomalies of its tidal response caused by horizontal non-uniformities. Finally, extremely small and stable drift of SG helps to study for non-tidal gravity variations caused by polar motion, vertical crustal movements etc.
The majority of the phenomena under consideration are spatially coherent, therefore its investigation is possible only on the global basis. The overall amount of installed (or planned) SGs is about 30. All the devices are working in the framework of international GGP project (Global Geodynamical Project). During the last years new SGs were installed in Chile, Chechien, Sweden, South Korea, India and Taiwan.

Top

Modelling of reactive fluid transport in deformable porous rocks
V.M. Yarushina and Y.Y. Podladchikov University of Oslo, Physics of Geological Processes, Oslo, Norway (viktory@fys.uio.no)

One outstanding challenge in geology today is the formulation of an understanding of the interaction between rocks and fluids. Advances in such knowledge are important for a broad range of geologic settings including partial melting and subsequent migration and emplacement of a melt into upper levels of the crust, or fluid flow during regional metamorphism and metasomatism. Rock-fluid interaction involves heat and mass transfer, deformation, hydrodynamic flow, and chemical reactions, thereby necessitating its consideration as a complex process coupling several simultaneous mechanisms. Deformation, chemical reactions, and fluid flow are coupled processes. Each affects the others. Special effort is required for accurate modelling of the porosity field through time. Mechanical compaction of porous rocks is usually treated under isothermal or isoentropic simplifying assumptions. However, joint consideration of both mechanical compaction and reactive porosity alteration requires somewhat greater than usual care about thermodynamic consistency. Here we consider the modelling of multi-component, multi-phase systems, which is fundamental to the study of fluid-rock interaction. Based on the conservation laws for mass, momentum, and energy in the form adopted in the theory of mixtures, we derive a thermodynamically admissible closed system of equations describing the coupling of heat and mass transfer, chemical reactions, and fluid flow in a deformable solid matrix. Geological environments where reactive transport is important are located at different depths and accordingly have different rheologies. In the near surface, elastic or elastoplastic properties would dominate, whereas viscoplasticity would have a profound effect deeper in the lithosphere. Poorly understood rheologies of heterogeneous porous rocks are derived from well understood processes (i.e., elasticity, viscosity, plastic flow, fracturing, and their combinations) on the microscale by considering a representative volume element and subsequent averaging of microscopic constitutive laws. Micromechanical and thermodynamic modelling is performed in such a way that the consistency of the obtained rheology and thermodynamically admissible closed system of equations with the exact Gassman’s relationship and Terzaghi effective stress law in the simplified case of poroelasticity is guaranteed.

In such environments as subduction zones or mid-ocean ridge, metamorphic rocks exhibit a lack of chemical homogenisation. Geochemistry suggests that in order to produce chemical heterogeneity, the fluids generated during high-pressure metamorphism must have been strongly channelled. The following three major mechanisms of fluid flow focusing have been proposed: fluid flow in open fractures and two different types of flow instabilities that do not require the pre-existing fracture network. Of the latter, the first represents a purely mechanical instability of Darcian flow through the deformable porous rock while the second is reactive infiltration instability. Obtained model is applied to a study of mechanical instability of Darcian flow through the porous rock (porosity wave) that can nucleate from perturbations to an initially uniform porosity and grow by draining fluid from the background porosity. The influence of different rock rheologies on the properties of porosity waves is investigated.

Top

The Fate of the Slabs Interacting with a Viscosity Hill in Mid-Mantle.
David A. Yuen, University of Minnesota, U.S.A.
Gabriele Morra, University of Sydney, Australia
Larry Boschi and Paul J. Tackley, Swiss Federal Institute of Technology, Switzerland

In the last two decades it has been proposed several times already that a viscosity hill might be present in the middle of the mantle, at a depth roughly between 1200 km and 2000 km. Such a viscosity obstacle would have a significant geodynamical consequence . Many tomographical models display strong signals of the presence of “fast” material lying at mid mantle depths and a spectral analysis for seismic tomography displays a very clear transition for degree up to around 16 at a 1500km depth. Furthermore recent works, both theoretical and experimental, on the high to low spin transition for periclase, have suggested that the high-spin to low-spin transition of Fe++ might lie at the heart of all these observations, because of the radical change in the electronic environment in the 3d orbitals of Fe++ during this second-order transition .In order to verify the geodynamic plausibility of such hypothesis, we employ here a recently developed Multipole-accelerated Boundary Element (FMM-BEM) numerical approach for solving the momentum equation in the global spherical earth is employed for simulating the interaction of one or more sinking slabs with a mid-mantle density and viscosity smooth discontinuity. We have focussed on the complexities induced to the behavior of very large plates O(2000km-10000km), characteristic of the Kula-Farallon plate system 60 Myrs ago.. We drmonstrate that the combination of sphericity and a mid mantle density-viscosity discontinuity trigger an intricate set of slab morphologies. Lack of space at depth forces the subducting plate either to thicken or to fold. Furthermore, whether or not the plate can penetrate into the lowest mantle depending on the radial profile of density and viscosity, within a realistic range (viscosity up to 10 or 100 times more viscous of the rest of the mantle, differential densities in the range 0 to 4% toward the bottom of the lower mantle). We also show how the decrease of thermal expansivity would cause a considerable slowing of the sinking rate of the slabs in the lower mantle. We conclude that the high-spin to low-spin transition of Fe++ is indeed able to trigger a very wide spectrum of behaviors in the lower mantle and we propose a scenario through which in the long term this effect might create a set of stalling slabs, which is characterized by chemically depleted slabs in the mid mantle. This would chemically stratify partially the mantle at around these aforementioned depths.

Top

On hypothesis of hydrogen in the Martian core
V. N. Zharkov Schmidt Institute of Physics of the Earth, Russian Academy of Sciences, (Zharkov @ ifz. ru / Fax: +7-499-2556040)

The idea of the origin of the terrestrial planets from planetesimals oxidized to varying degrees was discussed by Dreibus and Waenke (1989). The composition of the Earth and Mars is considered to be a certain mixture of component A and B.

Component A. This substance is highly reduced. Iron and all siderophile elements are in the metal state, and even silicon is partially metallic. Protobodies consisting of the A component filled the feeding zone of the forming Earth.

Component B. This substance is highly oxidized and contains all elements, including volatiles, with abundances like those in meteorites of class C1. Iron and all siderophile and lithophile elements are present mainly as oxides. Component B constituted the protobodies in the zone of the contemporary asteroidal belt.

Dreibus and Waenke (1989) concluded from their estimates of the bulk composition of Mars that component A and B are mixed in this planet in the ratio 60:40, whereas the ratio 85:15 was reported for the Earth. In was inferred that the accretion of Mars was almost homogeneous, in contrast to the chemically inhomogeneous accretion of the Earth.

At present, a commonly adopted opinion is that Jupiter was the first among the planets to be formed (Zharkov, 1993). Owing to a strong gravitational field, young Jupiter strewed the remaining protobodies from its feeding zone. These protobodies, as well as resonant interactions, destroyed the feeding zone of the planet, which could be formed in the asteroidal belt, and substantially reduced the number of protobodies in the feeding zone of growing Mars, thus slowing down its growth. As a result, Mars acquired a mass an order of magnitude smaller than could be attained without the influence of Jupiter. It is Jupiter’s effect that is responsible for the two- component model of the formation of terrestrial planets.

In the planetary interior, the conditions for the dissolution of hydrogen in iron arise under high pressures and temperatures. The source of hydrogen is the reaction

Fe + H₂O → FeO + H₂.
As gas pressure grows, hydrogen starts dissolving in iron:

Fe + x/2H₂→ FeH_x.
The pressure at the center of Mars, according to the available models, is less than 450 kbar. Hydrogen dissolves in iron at relatively low pressures. Particularly, the molecular ratio H/Fe 0.2- 0.4 is attained near the melting point at pressures as small as several tens of kilobars, and the ratio H/Fe 1 is attained at moderate temperatures and pressures of ~70 kbar. The dissolution of hydrogen in iron reduces the density and appreciably lowers the melting point. We summarize the experimental data on the system Fe-S-H and advance arguments in favor of the presence of hydrogen in the core, estimate the influence of hydrogen on the physical parameters of the core.

The oxidized component B, mentioned at the beginning of the abstract, is assumed to have a composition equivalent to that of carbonaceous chondrites of class C1. The water content attains in such bodies ~7.3 wt %. Ahrens et al. (1989) studied in laboratory conditions the question as to what are the shock pressures at which the minerals containing volatiles (CO , H O, and SO ) begin to lose them. They found that consolidated minerals start losing volatiles at shock pressures within the range 30-50 GPa. Such shock pressures are achievable when the impactor has a velocity of ~2-3 km/s, colliding with a target of the same material. Consequently, the dehydration of planetesimals that belong to component B begins when the radius of growing proto-Mars is and attains 75% at (R being the radius of Mars).

To estimate the mass of H that can be buried in the interior of growing proto-Mars, we put r ≈ 0.5R and ρ ≈3.5g/cm³. Then one easily obtains ~2.4 X 10²³g of hydrogen. The mass of the iron-nickel (of Fe₀₉Ni_0.1) core comprises in model DW ~1.2 X 10²⁶g. Therefore, if all the buried hydrogen ultimately becomes a constituent of the Martian core, we obtain the composition (Fe₀₉Ni_0.1)₀₉H_0.1, where the subscripts mean molecular fractions. This estimate can be considered the lower bound for the hydrogen content in the Martian core. Under the assumption that component B has the same composition as C1 chondrites (~7.3 wt % of H₂O), we produced in the process of the accumulation of Mars is ~2X10²⁴g of H₂. Thus, the upper bound for the hydrogen content in the iron core of Mars, not achievable in practice, is

The estimates given above were obtained with the data on water content of about 7.3 wt. % in C1 chondrites from the paper by Dreibus and Waenke (1989). If one takes into account recent estimates of water countent (Lodders and Fegley, 1998) the estimates of hydrogen content in the Martian core should be increased by a factor of 2.5-3. Then, the estimates of hydrogen in the core as H_0.5 and H_0.7 look quite reasonable. Zharkov(1996) found that addition of the molecular concentration x=0.1 of hydrogen to the iron core of Mars reduces its density by about 0.16 g/cm³

References
Ahrens, T.J., O’Keefe, J.D., Lange, M.A., 1989. Formation of atmospheres during accretion of the terrestrial planets. In: Attreya, S.K., Pollack, J.B., Matthews M.S.(Eds.). Origin and Evolution of Planetary and Satellite Atmospheres. Univ. Arizona Press, pp.328-385.
Dreibus G., Waenke, H., 1989. Supply and loss of volatile constituents during the accretion of terrestrial planels. L.S..,pp.268-288.
Lodders, K., Fegley, B., 1998. The Planetary Scientist’s Companion. Oxford University Press.
Zharkov, V.N.,1996. The internal structure of Mars: a key to understanding the origin of terrestrial planets. Solar Syst. Res. 30,456-465.
Zharkov, V.N., 1993. The role Jupiter in the formation of planets. Geophys. Monogr. 74, IUGG,
Am.Geophys. Union 14, 7-17.

Top

Dynamics and Implications of 3-D hydrous thermal-chemical plumes in oceanic subduction zones
Guizhi Zhu¹, Taras Gerya^1,2, Dave A.Yuen³, James A.D. Connolly¹
1. Institute of Geophysics, Swiss Federal Institute of Technology (ETH Zurich), CH-8092 Zurich, Switzerland
2. Adjunct Professor of Geology Department, Moscow State University, 119899 Moscow, Russia
3. University of Minnesota, Minneapolis, MN 55455-0219,USA

Magmatism and seismic data in subduction zones show that subduction is inherently three-dimensional (3-D). This 3-D character has traditionally been explored with thermo-mechanical models, i.e. only with thermal buoyancy. Here we expand on these studies by incorporating the thermal-chemical effects of dehydration/hydration processes and water transport above slabs in a 3-D petrological (i.e. chemical-) thermo-mechanical model for intra-oceanic subduction. The simulations were carried out with I3ELVIS code (2009), which is based on a multigrid conservative finite-difference approach combined with marker-in-cell methods. Our numerical simulations show that three types of plumes occur above the slab-mantle interface: 1) finger-like plumes that grow atop sheet-like trench-parallel structure; 2) Ridge-like structures perpendicular to the trench; 3) periodic wave-like structures propagating upwards along the upper surface of the slab that form zig-zag patterns parallel to the trench.

The viscosity of the plume material is the main factor controlling the development of the plumes. Our results show that lower viscosity of the partially molten rocks facilitates the Rayleigh-Taylor like instabilities at small wavelengths.

Colored streamlines in the mantle wedge reveal the small-scale circulations associated with these small-scale instabilities (e.g., see Fig.1a and 1b). It shows the importance of thermal-chemical buoyancy in 3-D mantle wedge. These small-scale motions obfuscate long wavelength fluid motions, as advocated by the trench-parallel flow model to explain complex seismic anisotropy in subduction zone proposed by e.g. Kneller (2005) and Silver ( 2008). Recent findings by Faccenda et al. (2008) and Healy et al. (2009) also would take exceptions to this trench-parallel flow orthodoxy.

Fig.1. Dynamics of development of small-scale flow field (stream lines) during the growth of the shell-like structure (wave). Stream lines are obtained by subtracting the average long-term slab retreat velocity of 2 cm/yr. Colors of stream lines indicate the magnitude of flow velocity: red=large, blue=small. Note splitting of convection pattern to several vortexes visible below the wave at the late stage of the process. This splitting is presumably related to the growth of trench-orthogonal ridges behind the wave. We emphasize here that these small-scale flows would frustrate the formation of any long-wavelength trench-parallel flow proposed by Silver (2008)

Reference: 1. Gerya, T. V. (2009), Introduction to numerical geodynamic modelling, Cambridge University Press (in press), Chapter 14 and 15.
2. Kneller, E.A., van Keken, P.E., Karato, S., et al. (2005), B-type olivine fabric in the mantle wedge insight from high-resolution non-Newtonian subduction zone models. Earth and Planetary Science letters, 237,781-797.
3. Long, A. D. and Silver, P.G. (2008), The subduction zone flow field from seismic anistotropy : A global view, Science, 319,315-318.
4. Faccenda, M, Burlini, L, Gerya, T.V., et al. (2008), Fault-induced seismic anisotropy by hydration in subducting oceanic plate, Nature, 455,1097-1100.
5. Healy, D., Reddy, S. M., Timms, N., et al. (2009), Trench-parallel fast axes of seismic anisotropy due to fluid-filled cracks in subducting slabs. Earth and Planetary Science letters, 283,75-86.

Top

Crustal and deep mantle structure of North Atlantic MOR
Kalinina A. V., A. A. Baranov
Schmidt Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russia (kalinina_av@mail.ru, baranov@ifz.ru)

Many people in different parts of earth sciences use now global CRUST2.0 (Bassin et al., 2002) model with resolution 2x2 degree. But for present days this model stays rather poor and has not enough resolution for accurate calculations. We begin to build complex geophysical model for North Atlantic MOR using different geophysical data. It was used data of deep seismic reflection, refraction, receiver functions and gravity studies from published papers. The existing data were verified and crosschecked. North Atlantic MOR is rather good investigated and geophysical data about it`s building help to build plate tectonics theory. Using new geophysical data some conclusions can be done about possibility of existence some blocks of continental crust in oceanic basins and either in MOR areas.

Also seismic and gravity data shows us crustal density asymmetry of opposite sides of North Atlantic MOR. Seismic and gravity data shows us density of upper crustal layers is about 2.7- 2.8 g/sm³ on east side but on west side this is approximately about 2.3- 2.4g/sm³. (density was received using recalculation from Vp and Vs velocities and using gravity modeling data) (Budanov et al., 2000). Also geophysical investigations of some transform faults (Kane, Hayes, Vema and others) shows big lateral inhomogeneities in their deep structure and properties of the crust and upper mantle. When we understand this we decide to improve old crustal model using seismic data which stay available during last years.

References:
1. Bassin et al., The Current Limits of Resolution for Surface Wave
Tomography in North America, // EOS Trans AGU. 2000. 81(48), Fall Meet. Suppl., Abstract
F897, http://mahi.ucsd.edu/Gabi/rem.html

Top

Gravity stresses in the Central Asia crust.
V. Pogorelov, A. Baranov
Schmidt Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russian Federation (vpogorelov@list.ru, baranov@ifz.ru)

The Southern and Central Asia is a tectonically complex region which characterized by the great collision between the Asian and Indian plates. Its tectonic evolution is strongly related to the active subduction process along the Pacific border. Stress investigation in the continental crust is a very important problem not only for science but also for the practical purposes. There are four main factors which produce tectonic stresses: gravity anomalies of the crust, density inhomogeneities, deformation from area with intraplate collision, residual elastic deformationsand underthrust stresses conditions from convective mantle. We present the stress model of the crust and lithosphere for the Central and Southern Asia on the basis of the finite element modeling and use gravity anomalies and density inhomogeneities of the crust as loading mechanism.

Basing on new crustal model ASCRUST-08 with resolution 1º×1º, which was built using new seismic data [Baranov, 2008], it were analyzed three vertical density profiles: 84, 89 and 94 W.E, (fig. 1, fig. 2).

Figure.1. Relief as a base map and Moho map (contours) using model (ASCRUST-08) for Central Asia.

Figure.2. Mass distribution on Moho border.

For the crust we take the elasto-plastic rheology with Drucker-Prager criterion. In the lithosphere the elasto-plastic model with von Mises criterion is assumed. We investigated stresses which are produced by the crustal density inhomogeneities and surface relief. The calculations are done using the U-WAY finite element code developed at the Institute of Applied Mechanics Russian Academy of Sciences. (similar to the Nastran program) Density inhomogeneities are based on the AsCRUST-08 crustal model (Baranov, 2008), which has resolution of 1 x 1 degree. AsCRUST-08 was built using the data of deep seismic reflection, refraction and receiver functions studies from published papers. The complex 3D crustal model consists of three layers: upper, middle, and lower crust. Besides depth of the boundaries, we provided average P-wave velocities in the upper, middle and lower parts of the crystalline crust and sediments. The seismic P-velocity data was also recalculated to the densities and the elastic modules of the crustal layers using the rheological properties and geological constraints. (Barton, 1986; Pauselli et al., 2003.)

Strength parameters of rocks strongly depend on temperature, tectonic and fluid pressure. Fluid pressure can reduce resistance forces in faulting rock, tectonic pressure increases these forces.

References
1. Baranov A. A. Complex 3D crustal model of Central and Southern Asia //Electronic journal «Vestnik RAS» №1(26)2008 (URL: http://www.scgis.ru/russian/cp1251/h_dgggms/1-2008/informbul-3_2008/cw-3.pdf (in Russian)

2. Barton, P.J., 1986. The relationship between seismic velocity and density in the continental crust—a useful constrain. Geophys. J. R. Astron. Soc. 87, 195– 208.

3. Pauselli C. and C. Federico «Elastic modeling of the Alto Tiberina normal fault (Central Italy): geometry and litological stratification influences on the local stress field». Tectonophysics, No. 347, 2003 –p. 99-113.

Top