Бесплатный автореферат и диссертация по наукам о земле на тему

Комплексное компьютерное моделирование геохимических объектов на примерах двумерных моделей коллизии плит, магматогенно-гидротермальной системы и зоны субдукции

ВАК РФ 25.00.09, Геохимия, геохимические методы поисков полезных ископаемых

Автореферат диссертации по теме "Комплексное компьютерное моделирование геохимических объектов на примерах двумерных моделей коллизии плит, магматогенно-гидротермальной системы и зоны субдукции"

■П с

На правах рукописи

ВАСИЛЬЕВ Владимир Игоревич

КОМПЛЕКСНОЕ КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ГЕОХИМИЧЕСКИХ ОБЪЕКТОВ НА ПРИМЕРАХ ДВУМЕРНЫХ МОДЕЛЕЙ КОЛЛИЗИИ ПЛИТ, МАГМАТОГЕННО-ГИДРОТЕРМАЛЬНОЙ СИСТЕМЫ И ЗОНЫ СУБДУКЦИИ

Специальность 25.00.09 - геохимия, геохимические методы поисков полезных ископаемых

Автореферат

диссертации на соискание ученой степени кандидата геолого-минералогических наук

ИРКУТСК 2009

003481987

Работа выполнена в Геологическом институте СО РАН (ГИНСОРАН)

Научные руководители:

доктор геолого-минералогических наук

Жатнуев Николай Сергеевич доктор геолого-минералогических наук

Чудненко Константин Вадимович

Официальные оппоненты: доктор химических наук

Таусон Владимир Львович доктор геолого-минералогических наук

Скворцов Валерий Александрович

Ведущая организация: Институт Земной коры СО РАН

Защита состоится ноября 2009 г. в часов на заседании диссертационного совета Д 003.059.01 при Институте геохимии им. А.П. Виноградова СО РАН по адресу: 664033, г. Иркутск, ул. Фаворского 1а.

С диссертацией можно ознакомиться в научной библиотеке Института геохимии им. А.П. Виноградова СО РАН, по адресу: 664033 Иркутск, ул. Фаворского 1а.

Автореферат разослан

о10 октября 2009 г.

Ученый секретарь диссертационного совета

к.г.-м.н.

Г.П. Королева

ВВЕДЕНИЕ

Актуальность работы. Прогресс геохимии зависит в первую очередь не от дальнейшего расширения банка аналитических данных, а от синтеза всей совокупности полевой, аналитической и экспериментальной информации, что наиболее перспективно при использовании компьютерного моделирования. Развитие методов термодинамического компьютерного моделирования в геохимии состоит в преодолении главного ограничения аппарата химической термодинамики - отсутствия пространственно-временных координат (Шарапов, 1986, 2007). В работе представлен подход к решению этой задачи, а также затрагиваются кардинальные проблемы моделирования на ЭВМ в геохимии (Карпов, 1981).

Цель исследования состоит в разработке методики комплексного компьютерного моделирования геохимических систем на основе объектно-ориентированного подхода и апробации ее на тематических моделях геохимических объектов. Задачи работы были поставлены в соответствии с методологическим характером выбранной цели и сформулированы следующим образом:

1) анализ и модернизация существующих методов и средств физико-химического моделирования;

2) наработка опыта в построении разнотипных компьютерных моделей и разработке прикладных программ, необходимых при их построении;

3) разработка методики моделирования геохимических объектов и апробация ее на примерах модели-схемы коллизии плит, модели магматогенно-гидротермальной системы и комплексной модели зоны субдукции;

4) анализ, интерпретация и оценка результатов моделирования геохимических систем.

Защищаемые положения:

1. Предложенная методика комплексного компьютерного моделирования позволяет в рамках универсального объектно-ориентированного подхода строить модели самых различных геохимических объектов и оценивать их по предложенным критериям корректности: описательному, функциональному и событийному.

2. Модель-схема коллизии плит, разработанная по предложенной методике и реализованная с помощью прикладного программного продукта Vladi Collision, позволяет, варьируя начальными и граничными условиями, рассчитывать и визуализировать различные варианты геометрической эволюции модельной системы с учетом распределения температурных полей разреза.

3. Модель магматогенно-гидротермалъной системы отражает основные закономерности массопереноса и минералообразования в гидротермальных системах областей развития современного вулканизма.

4. Разработанная комплексная компьютерная модель зоны субдукции, представленной в виде двумерной метасистемы за счет разбиения на более чем 4000 подсистем, подтверждает гипотетические положения о Р-Т-условиях нахождения свободного флюида и его эволюции при подъеме от дегидратирующейся плиты к поверхности.

Научная новизна:

1) Впервые в рамках методики единого подхода сведены вместе четыре аспекта моделирования геохимических объектов: геометрический, физический, физико-химический и динамический, что позволило увеличить параметры корректности при переходе от эмпирических моделей к компьютерным.

2) Доказана принципиальная возможность использования при моделировании геохимических объектов объектно-ориентированного подхода, позволяющего максимально наглядно и удобно представлять объекты в терминах конструктивных схем единой концепции.

3) Получена полная двумерная физико-химическая и геохимическая картина для модели зоны субдукции: рассчитаны более четырех тысяч равновесных минеральных парагенезисов с раствором и газовой фазой, температуры и давления в каждой точке модельной плоскости. Определены области дегидратации и плавления погружающейся плиты, рассчитан перенос флюидом более 90 зависимых компонентов с глубин -100 км к поверхности.

4) Построена физико-химическая модель на материале Северо-Парамуширской магматогенно-гидротермальной системы. Определены Р-Т-параметры, рассчитаны равновесные химические составы на пути конвективного движения флюида, выявлены формы переноса более 150 компонентов водного раствора и характеристики изменения вмещающих пород.

Фактическим материалом для данной работы явились общепризнанные литературные данные по химическому и минеральному составу горных пород, их физическим, петрофизическим, термодинамическим свойствам, представленные в различных справочных, монографических и журнальных изданиях. В работе также использовался фактический материал, полученный автором в полевых работах сезонов 1995-2008 гг., в процессе научно-исследовательских работ по теме «Условия и механизмы формирования благороднометального оруденения (Аи, ЭПГ) традиционных и новых типов в Саяно-Байкало-Муйском сегменте Цен-

трально-Азиатского складчатого пояса (по природным и экспериментальным данным)» (№ гос. per. 01.2.007 05169) и при финансовой поддержке Президиумов СО и ДВО РАН (проект № 117 - 09-II-CO-08-006). Кроме этого, фактические данные по химическому составу, петрофизи-ке и теплофизике учитываемых в моделировании природных сред были любезно предоставлены непосредственно Н.С. Жатнуевым, В.И. Гуни-ным (ГИН СО РАН), И.К. Карповым, К.В. Чудненко, В.А. Бычинским (ИГХ СО РАН), С.Н. Рычаговым (ИВиС ДВО РАН) и другими учеными. Практическая ценность работы. Представлены методологически единые модели таких категориально различных объектов, как процесс коллизии плит и объект термодинамической системы в толще породы. Разработанные модели для Северо-Парамуширской гидротермально-магматической системы и зоны субдукции показывают практическую пригодность методики для получения как качественных, так и количественных параметров геохимических объектов. Предлагаемые критерии корректности рекомендуется использовать для критической оценки и сравнения компьютерных моделей при решении задач по эволюции вещества в геохимических процессах. Основные положения работы использованы при разработке и чтении курсов лекций, проведении практических занятий и полевых практик на Кафедре геологии Бурятского государственного университета в 2000-2009 гг.: «Информатика», «Компьютерное моделирование», «Геодезия и картография», «ЭВТ в геологии», «Основы физико-химического моделирования в геохимии». Публикации и апробация работы. Основные положения диссертации опубликованы в 18 печатных работах, а также докладывались на Школе-семинаре российских делегатов XXXI Международного Геологического Конгресса (НИС «Академик Иоффе», Атлантика, 2000), XIX Всероссийской конференции «Строение литосферы и геодинамика» (Иркутск, 2001), III Всероссийском симпозиуме с международным участием «Золото Сибири и Дальнего Востока» (Улан-Удэ, 2004), III Всероссийском симпозиуме по вулканологии и палеовулканологии (Улан-Удэ,

2006), Всероссийской конференции с иностранным участием в честь 50-летия СО РАН и 80-летия чл.-корр. РАН Ф.П. Кренделева (Улан-Удэ,

2007), I международной конференции «Граниты и эволюция Земли» (Улан-Удэ, 2008), IV Всероссийском симпозиуме по вулканологии и палеовулканологии «Вулканизм и геодинамика» (Петропавловск-Камчатский, 2009), Научных совещаниях «Геодинамическая эволюция литосферы Центрально-Азиатского подвижного пояса (от океана к континенту)» (Иркутск, 2008, 2009), ежегодных Научных сессиях Геологического института СО РАН (1998-2009) и др. Различные варианты авторских программных продуктов переданы и используются в лаборато-

риях Геологического института СО РАН (г. Улан-Удэ) и Института геохимии СО РАН (г. Иркутск).

Объем работы. Диссертация состоит из введения, пяти глав, заключения, списка литературы и приложения общим объемом 120 страниц печатного текста, 16 таблиц и 36 иллюстраций. Список использованной литературы включает 159 наименований.

Благодарности. Автор считает своим долгом выразить сердечную благодарность в первую очередь своим научным руководителям - докторам геолого-минералогических наук Игорю Константиновичу Карпову, Николаю Сергеевичу Жатнуеву и Константину Вадимовичу Чудненко за терпеливую плодотворную помощь на протяжении всех лет работы над диссертацией. Автор душевно благодарит своих учителей - д.г.-м.н. Дмитрия Ивановича Царева, д.г.-м.н. Анатолия Георгиевича Миронова, к.г.-м.н. Валерия Алексеевича Бычинского, Владимира Ивановича Гу-нина и всех тех, кто на протяжении многих лет помогал вначале осваивать, а затем развивать и внедрять компьютерное моделирование в практику геологических исследований.

Глава 1. Состояние проблемы

Геометрический аспект. Степень геометрической детализации модели основывается на применении метода декомпозиции. Объект при этом представляется как система агрегатов, которые, в свою очередь, состоят из узлов, а те - из деталей (Быков, 2002). Замена объемной задачи на ряд плоских долго оставалась единственным способом решения и обычно приводит к приемлемым результатам и сегодня.

Физический аспект. Численные методы исследования математической модели предполагают создание компьютерной программы, моделирующей изучаемое явление. Использование компьютера - часто единственный способ получения следствий го математической модели. Компьютерный физический эксперимент фактически является физическим экспериментом над математической моделью исследуемого объекта, проводимого с помощью компьютера (Майер, 2009).

К классическим схемам данного аспекта моделирования могут быть отнесены: модель охлаждения интрузива JI.M. Кэтлса, модель теп-лопереноса A.B. Кирюхина и В.М. Сугробова, комбинированные модели Д.В. Гричука, М.В. Борисова и A.B. Тутубалина и др.

Физико-химический аспект основывается на применении методов термодинамического моделирования. В настоящее время существует два подхода к решению проблемы расчета химического равновесия (полного, промежуточного, метастабильного, частичного) в сложных многокомпонентных, многофазовых и многоагрегатных природных системах. Первый основан на формализме стехиометрических уравнений

реакций и константах их равновесия - расчет по реакциям, второй - на постановке и решении задач химического равновесия как задач выпуклого программирования - метод минимизации. Расчет по реакциям использовался химиками и технологами в докомпьютерную эру развития науки. Среди зарубежных геохимиков программы, в основу которых положен метод констант равновесия, получили широкое распространение и применение в изучении процессов взаимодействия вода - горные породы главным образом благодаря пионерским работам Г. Хельгесона (Не^евоп, 1967, 1968, 1969, 1970; Не^евоп & а!., 1969, 1970). Им к началу 1970-х годов были разработаны математическая модель, вычислительный алгоритм, методы формирования базы термодинамических данных с участием компонентов водных растворов электролита; создана компьютерная программа; а, главное, показан принципиально новый подход к моделированию физико-химических процессов с учетом их необратимости на примере модели образования метасоматической зональности гидротермальных месторождений.

Постановка задач химического равновесия в формулировке выпуклого программирования имеет определенные преимущества: в каждом варианте решения оно дает на выходе, по крайней мере, в два раза больше базовой термодинамической информации. Помимо мольных количеств независимых компонентов (или других единиц содержания), включая фазовый состав, мы, используя минимизацию, получаем численные значения химических потенциалов независимых компонентов и детальную характеристику решений, с помощью которых молено сделать вывод о достижимости глобального минимума и внутренней согласованности исходных термодинамических данных тех веществ, которые вошли в оптимальное решение. Имеется и ряд других возможностей термодинамического моделирования в формулировке выпуклого программирования, недоступных методу реакций, что объективно показано во многих работах (Карпов, 1981; Карпов и др., 1981, 1987, 1991; Бакшеев, Карпов, 1984, 1988; Чудненко и др., 1987, 1988, 1999; КиНк, 2000,2002).

Глава 2. Методика моделирования

Компьютерные модели реальных геохимических объектов должны отображать основные закономерности поведения исследуемого природного прототипа на основании применяемых вычислительных методов и алгоритмов с учетом как можно большего количества наблюдаемых, экспериментально определенных или предполагаемых эмпирических характеристик. Будем понимать под термином объект - совокупность характеристик рассматриваемого явления. Каждый объект имеет следующие характеристики:

свойства - размеры, состав, положение в пространстве и т.д.; методы - способы изменения свойств объекта; события - внешние воздействия, на которые объект может реагировать определенными для него методами.

Следовательно, построение модели объекта - это описание необходимых свойств, потенциальных событий и методов, которыми объект будет реагировать на эти события. Точность и полнота модели определяется отношением учтенных характеристик к их общему количеству и может быть оценена тремя параметрами, характеризующими ее корректность: описательным, событийным и функциональным.

Допустим, что физический (природный) объект обладает

множеством свойств Р о {1...р} и множеством событий Е з {1___е}, на

которые он может реагировать множеством методов Мэ {1 ...т}. Модель этого объекта может учитывать не более рт свойств, ет событий и тт методов, где, как правило, рт<<р) ет«е и тга«т. Тогда описательная корректность модели КР, равная отношению рт/р учтенных в модели свойств к их реальному количеству, будет определять статическое соответствие модели реальному объекту. Событийная корректность КЕ = ет/е оценивает в модели учет влияния внешних факторов, а функциональная корректность Км = Шп/т позволяет оценить соответствие динамики модели поведению реального объекта (эволюции объекта). Стремление параметров корректности к единице свидетельствует о том, что сложность модели приближается к сложности природного объекта.

Предлагаемая комплексная методика моделирования позволяет максимально увеличить величины параметров корректности на этапе перехода от эмпирической модели к компьютерной способом поэлементного учета характеристик. Она включает четыре взаимосвязанных этапа: геометрический, физический, физико-химический и динамический.

Геометрический этап. Начиная создание модели, необходимо построить геометрическую схему, которая должна наглядно передавать свойства объекта, как в пространственном, так и в понятийном аспекте, обеспечивая высокий уровень описательной корректности. Основные шаги этапа заключаются в следующем. Во-первых, необходимо выбрать геометрическую размерность модели. Во-вторых, любая достаточно сложная модель представляет собой мегасистему, состоящую из конечного количества подсистем, геометрические соотношения между которыми также необходимо предусмотреть. В общем случае это производится разбиением моделируемого пространствана геометрические фигуры, каждая из которых соответствует определенной подсистеме со своими физическими и химическими параметрами. В-третьих, для со-

блюдения соответствия модель-объект на макроуровне, необходимо произвести дискретизацию пространства, разделив все имеющееся модельное пространство на зоны, состоящие из конечного количества подсистем. Данные зоны должны соответствовать принципиальным природным обособлениям в составе моделируемого объекта.

Физический этап предусматривает расчет физических параметров каждой подсистемы, а также физические аспекты взаимодействия подсистем между собой. Рассмотрим суть физического этапа на примере моделирования тепломассопотоков в гидротермальной системе. Набор физических характеристик подсистем в этом случае должен включать массу, плотность, объем, пористость, а также теплофизиче-ские свойства их сред. Начальными и граничными условиями физического этапа являются начальные градиенты температуры и давления, проницаемость границ и режим теплообмена подсистем и модели в целом и т.д. В качестве расчетных параметров, как правило, выступают температура и давление.

Предлагается следующий метод расчета распределения температурных полей в двумерном модельном пространстве. Сначала рассчитываем кондуктивную составляющую теплопереноса на основе стандартных дифференциальных уравнений (Аверкин, Шарапов, 1986). Коэффициент теплопроводности при температуре Т для пород магматического состава можно рассчитывать по зависимости Тихомирова (Теп-лофизические свойства..., 1987). Для расчета теплоемкостей подсистем удобно использовать стандартное ступенчатое уравнение, описывающее зависимость теплоемкости от температуры Т (Чудненко, 2005). Также должна учитываться поправка на изменение теплоемкости, обусловленная эффектом разупорядочения в зависимости от температуры для подсистем, содержащих ряд минералов, таких как калиевый полевой шпат, доломит, геленит (Чудненко, 2005). При необходимости после расчета стационарного кондуктивного распределения температуры можно перейти к расчету конвективного теплопереноса, имеющего место при движении поровых растворов. Для пористой водонасыщенной среды в двумерной модели система дифференциальных уравнений определена в работах (Тихонов, Самарский, 1975; Жатнуев, Гунин, 2000) при следующих допущениях: поток через трещиноватую или пористую формацию подчиняется закону Дарси; движение раствора - двумерное, что дает возможность учета анизотропии среды по проницаемости; раствор и вмещающие породы находятся в состоянии термического равновесия. Для решения системы уравнений, описывающих процессы конвективного теплопереноса в двумерной модели, необходимо задать соответствующие начальные и граничные условия и использовать численные ме-

тоды, позволяющие рассчитать значения искомых функций в каждой подсистеме модельной плоскости. В нашем случае для численной реализации поставленной задачи можно использовать абсолютно устойчивый локально-одномерный метод решения (Тихонов, Самарский, 1975) с последующим методом установления с нижней релаксацией (Тарунин, 1990).

Литостатическое давление в моделируемых подсистемах можно рассматривать как линейную функцию количества геометрически вышележащих подсистем и плотностей их сред. Гидростатическое давление раствора (флюида) входит в набор учитываемых физических величин в случае наличия в составе подсистемы протяженных по вертикали трещин, где раствор находится в объемной фазе, и рассчитывается аналогично.

Физико-химический этап предусматривает расчет исходных равновесных составов подсистем, которые рассматриваются как физико-химические резервуары, изначально не взаимодействующие друг с другом. Входными данными для таких расчетов будут являться рассчитанные на предыдущем этапе температура и давление в каждом резервуаре, а также набор и концентрации независимых компонентов. Список независимых компонентов формируется, исходя из общей тематики моделирования. За основу принимается набор основных петрогенных элементов, входящих в состав оксидов стандартного силикатного анализа: "Л, А1, Бе, Мп, М^, Са, 1Ма, К, Р, О, Н, с добавлением компонентов, поведение которых необходимо исследовать в модели. Концентрации независимых компонентов каждой зоны (вектор Ь) определяются следующим образом. Сначала выбирается соответствующий зоне эмпирический вещественный состав, полученный в результате анализов проб или по литературным данным. В случае моделирования глобальных объектов (например, в случае геодинамических моделей) необходимо выбирать наиболее общий исходный состав (статистически усредненный из как можно более представительной выборки). В случае же построения модели локального геохимического объекта с достаточно изученным вещественным составом, последний должен быть отражен наиболее полно. К базовому составу добавляются содержания необходимых в модели компонентов в той химической форме, в какой они анализировались.

Если температура и давление каждой подсистемы индивидуальны, то набор независимых компонентов должен быть постоянным для всей метасистемы в целом, тогда как их концентрациями можно варьировать, исходя из принятых химических составов выбранных зон.

Расчет равновесных составов резервуаров предлагается про-

изводить методом минимизации потенциала Гиббса - С(Т, Р) на основе программного обеспечения, содержащего достаточно полные базы термодинамических данных зависимых компонентов, необходимых в модели. В качестве программного средства, реализующего физико-химический аспект моделирования, мы остановили свой выбор на программном комплексе «Селектор», с помощью которого доступно решение широкого круга задач, связанных с изучением физико-химических процессов в геохимии, петрологии и других областях. Модуль формирования моделей мегасистем позволяет решать задачи, связанные с исследованиями физико-химической эволюции взаимодействующей совокупности резервуаров по координате времени (Чудненко, 2007). Результатом этого этапа моделирования должны являться расчетные равновесные химические составы каждого резервуара (подсистемы).

Динамический этап. Начальными и граничными условиями здесь являются элементы дискретизации временной области. Главное из рассматриваемых событий при исследовании динамики объекта - это временные изменения, поэтому в первую очередь важно предусмотреть методы объекта как функции вектора времени. При этом учитываемые методы можно разделить на две категории по масштабу приложения: межзонные методы контролируют изменение свойств зон модели; межрезервуариые методы описывают взаимодействие резервуаров и изменение их свойств.

Первая категория иллюстрируется моделью-схемой коллизии плит, описанной в главе 3, вторая - моделями гидротермальной системы и субдукции в главах 4 и 5.

Глава 3. Модель-схема коллизии плит

В основу идеи модели положим простое распределение поступающего с движением плиты вещества в виде трех зон (рис. 1): треугольной верхней (гора), треугольной нижней (погружаемая в мантию часть) и прямоугольной центральной (мощность плиты). Модельная система симметрична относительно оси с.

Если скорость движения плиты равна {/, то площадь А51 поступающего за период времени Д/ плитного вещества равна АБ = Аа-М, где Аа = и-А( и М - мощность плиты. Для соблюдения условий подобия введем несколько коэффициентов, контролирующих форму и размеры зон.

Распределение вещества между треугольниками 1 я 2 будет контролироваться коэффициентом баланса Кв, который представляет собой расход поступающего вещества на треугольник 1 в долях от единицы. Соотношение высоты Я и основания I треугольника 1 контролирует коэффициент К у, который пропорционален вязкости вещества пли-

ты. Приращение длины основания треугольника 1 за время At рассчитывается по уравнению: д/, - 2KVKB AS + HL ^ а приращение высоты тре-

H-L

угольника 1 за время At - по уравнению: д// _ 2Л.'ВД5,(1- Прираще-

L + AL

ние площади нижнего треугольника 2 за время At будет равно: (1-А^ДО-

Форма нижнего треугольника (соотношение глубины D и длины основания В) будет контролироваться как коэффициентом вязкости Ку, так и коэффициентом экструзии мантии КЕ, зависящим от плотностей мантии и плиты. Приращение длины основания нижнего треугольника А В и глубины погружения AD за время A t рассчитываются по уравнениям: дд... 2KrKE(l - K„)AS + DB и др = 2(1 - KB)AS(\-KvKf )

D-B В + АВ

По достижении определенной высоты Нтах треугольника 1 (горы), стадия орогенеза сменяется стадией денудации, и появляется треугольная зона 3, представляющая собой сносимый (осадочный) материал. На этом этапе увеличение площади верхнего треугольника (7 + 3) идет только за счет увеличения длины стороны F треугольника 3. Как видим, все вышеприведенные методы изменения свойств принятых зон, выраженные через уравнения, содержат функцию времени. Модель учитывает два события: течение времени и достижение «горой» высоты Нтах. Вследствие изменения форм и размеров зон в процессе эволюции модели будет меняться и количество подсистем, входящих в состав каждой зоны.

Следовательно, мы должны в каждый моделируемый момент времени возвращаться на геометрический этап или делать это через какой-то временной шаг (Васильев, Жатнуев, 2007). Данная модель-схема была реализована автором в программный продукт Vladi Collision на языке С++..

Глава 4. Модель магматогенно-

гидротермальной системы Модель имеет вид проточного физико-химического реактора с расчетным начальным равновесным магматическим флюидом, который

I с

Рис. 1. Геометрическая схема распределения вещества при исследовании динамики коллизии двух плит.

в процессе линейного движения просачивается через конечное число резервуаров, каждый из которых имеет собственный химический состав, температуру и давление. Под «флюидом» в данной работе подразумевается группа подвижных фаз «раствор + газ». Все резервуары характеризуются единым набором независимых компонентов, исходные содержания которых для каждого резервуара задаются при постановке задачи. Задачей моделирования является расчет равновесных концентраций зависимых компонентов в каждом резервуаре, а также расчет состава флюида на выходе из последнего резервуара.

Модель

может усложняться притоком внешнего вещества в некоторые резервуары - в этом случае сначала рассчитывается равновесие резервуара с привнесенным веществом, а затем - с веществом флюида (рис. 2). Равновесные концентрации зависимых компонентов рассчитываются методом минимизации свободной энергии Гиббса. Для модели необходимы следующие исходные данные: температура, давление и химический состав вещества в магматическом очаге, химические составы каждого резервуара и Р-Т условия трассы движения флюида, заданные порезерву-арно или функцией.

В качестве примера рассмотрена конкретная модель подъема магматического флюида на материале, полученном при изучении Севе-ро-Парамуширской гидротермально-магматической системы. За вещественную основу модели был принят химический состав пород в разрезе скважины ГТТ-3 по данным исследования шлама и керна (Рычагов и др., 2002). Из данного состава в качестве независимых компонентов модели выбраны следующие элементы: А1, Бе, М§, Са, К, Б, С, Н и О. Добавлены рудные компоненты (Аи, Ag, Си, РЬ, Zn, Бп) и С1. Для

Рис. 2. Схема массопе-реноса модели (проточный реактор). Обозначения: МО - магматический очаг, Р1; Р2, Рм -физико-химические резервуары; Фи - исходный флюид, равновесный с веществом очага в его Р-Т условиях; Фк -конечный флюид на выходе из последнего резервуара; Фп - привносимое в некоторый резервуар вещество (данный параметр опционален, в случае привноса вещества сначала рассчитывается равновесие резервуара с привнесенным веществом, а затем - с веществом флюида).

последних взяты их кларковые содержания (Справочник по геохимии..., 1990). Анализы по скважине даны для 22-х точек по глубине (от 64 до 2500 м), следовательно, были приняты 22 резервуара с нумерацией снизу вверх (в направлении движения флюида). Масса каждого резервуара принята равной 1 кг. Состав первичного магматического флюида был рассчитан из равновесия чистой воды со справочным средним составом андезита (Богатиков и др., 1987) при 900°С и 772 бар, что в данных условиях соответствует глубине 3000 м. В резервуары с 20 по 22 (глубины 372 - 64 м) дополнительно вводился метеорный раствор, состав которого был усреднен из справочных литературных данных.

Температуры резервуаров интерполировались полиномиальной интервальной регрессией из данных (Белоусов и др., 2002). Лито-статическое давление в резервуарах рассчитывалось исходя из принятой средней плотности пород 2,57 г/см3. Петрофизические свойства пород в разрезе скважины ГП-3 были предоставлены С.Н. Рычаговым.

Зависимые компоненты были отобраны из баз данных ПК «Селектор». Это 170 конденсированных фаз, 155 компонентов водного раствора и 9 компонентов газовой фазы. Модель была рассчитана для подъема 1 кг флюида. Расчет показал, что в резервуарах с 1 по 16 (глубины от 2500 до 1090 м) флюид находится только в газовой фазе, а начиная с 20 резервуара (глубины от 372 до 64 м) - только в жидкой. На глубинах от 800 до 550 метров (резервуары 17 - 19) обе фазы сосуществуют. Характерно, что первичный флюид массой 1 кг попадает в первый резервуар (глубина 2500 м), и, как видно из диаграммы (рис. 3), его масса в следующем же резервуаре (2340 м) уменьшается на порядок, порода насыщается газами, в основном кислородом и соединениями серы.

Расчет показал, что в резервуарах с 17 по 22 (глубины от 800 до 64 м) существуют ненулевые концентрации 110 зависимых компонентов. Компоненты породообразующих элементов по поведению можно условно разделить на три группы. Концентрация в растворе компонентов

Масса флюида, кг

250

1250

1500

1750

2250

2500

Рис. 3. Расчетное изменение массы флюида (раствор + газ) по глубине. Маркерами показано расположение модельных резервуаров.

первой группы относительно высока на глубинах 800-700 м, а к поверхности уменьшается практически до полного отсутствия. Это справедливо для соединений Ы, Бе в хлоридной форме, N3 и К. Из рудных сюда можно отнести соединения Аи, Си и Н§. Ко второй группе отнесем компоненты с противоположным поведением - их концентрации в растворе существенно возрастают к поверхности. Это соединения Са, а из рудных - А.%, РЬ и Хп. Третья группа характеризуется максимальными концентрациями на средних для раствора глубинах - от 600 до 200 м. Ниже и выше этого интервала их концентрации незначительны. Это А1, Ре в свободной форме и М§, а из рудных элементов - Бп.

Ненулевые равновесные минеральные концентрации в различных резервуарах характерны для 35 конденсированных фаз. Рассмотрим рудные компоненты. Золото, серебро и ртуть значительных концентраций конденсированных фаз не образуют. Медь присутствует в виде халькопирита (интервал 372-64 м, до 0,031 масс. %), тенорита (интервал 1310-1090 м, до 0,015 масс. %), куприта (глубина 800 м, 0,011 масс. %) и халькозина (глубина 701 м, 0,012 масс. %). Единичное ненулевое содержание самородной меди зафиксировано на глубине 550 м. Цинк присутствует в виде цинкита (интервал 1646-550 м, до 0,029 масс. %) и сфалерита (интервал 372-64 м, до 0,042 масс. %). Олово на глубинах ниже 1300 м равновесно в самородном виде (до 0,009 масс. %), выше этой отметки - в виде касситерита (до 0,010 масс. %).

В целом, расчет показал пригодность данной концепции для моделирования массопереноса и минералообразования в магматогенно-гидротермальных системах. Результаты достаточно хорошо соответствуют литературным и природным данным (Васильев, 2009).

Глава 5. Модель зоны субдукции

Модель зоны субдукции является наиболее проработанной в соответствии с предложенной методикой, поэтому в данной главе дается наиболее полное описание последовательности этапов моделирования, начиная с постановки задачи.

Несмотря на длительную историю изучения вулканизма островных дуг и активных континентальных окраин, геотектонических позиций и глубинного строения зон его проявления, многие аспекты генерации магм и формирования геохимических характеристик магматических пород все еще остаются в значительной мере неопределенными. Ключевыми являются два основных вопроса (Авдейко и др., 2001):

1) что является источником магм, то есть «что плавится и где?»;

2) что является причиной плавления - дополнительное тепло, снижение давления или привнос флюидов?

Одни авторы (Добрецов и др., 2001) полагают, что плавится

субдуцируемая океаническая кора, другие (Kobayashi, 1981 и др.), - что плавится как океаническая кора, так и основание мантийного клина. Третьи считают, что плавится только вещество мантийного клина над субдуцируемой пластиной под воздействием воды и других летучих компонентов, отделяющихся от погружающейся океанической коры (Kushiro, 1975; Авдейко, 1994). Следовательно, поставленную задачу можно разделить на две составляющих:

1) Расчет физико-химически и математически обоснованных количественных соотношений Р-Т условий и равновесных концентраций минеральных парагенезисов всей модельной области.

2) Определение Р-Т-условий нахождения свободного флюида и расчет физико-химической эволюции вещества при его подъеме от этих областей к поверхности.

В свете данной проблемы была разработана и рассчитана модель распределения вещества и тепла для обстановки зон субдукции, включающая расчет подъема свободного флюида от зоны частичного плавления (~ 100 км) к поверхности.

Геометрический анализ и входные составы зон модели. Рассматриваемая в модели область представляет собой двумерный вертикальный разрез вкрест простирания зоны субдукции на глубину 100 км и горизонтальной протяженностью 40 км. Модельное пространство было разбито более чем на 4000 подсистем размером 1x1 км, для каждой из которых рассчитаны температура, давление и химический состав. Целые количества подсистем были объединены в условно принятые зоны, каждая из которых характеризуется собственным химическим составом, плотностью и пористостью. Принципиальное расположение зон заимствовано из работ (Winter, 2001; Richards, 2003) и представлено на рис. 4.

Согласно положению методики о предпочтительности наиболее общих исходных данных, входные химические составы принятых зон модели определялись на основании пропорционального соединения в них стандартных составов оболочек (слоев) континентальной, субконтинентальной, океанической коры и астеносферы. Такой «концептуальный состав» модельных зон приведен в таблице 1. Рассмотрим выделенные зоны.

Зона A («asthenosphere») имитирует верхнюю часть мантии валового пиролитового состава. Состав зоны SML («sub-continental mantle lithosphere») усреднен из литературных данных по «континентальному базальтовому» слою. Зона СС («continental crust») соответствует «гранитному» слою континентов. Принято, что континентальные зоны SML и СС не содержат метеорных вод ниже глубины 100-150 метров.

то есть, применительно к условиям модели, не содержат воды в объемной фазе. Зоны ОС («oceanic crust») и S («sediments») представляют соответственно океаническую кору и аккреционный клин морских и вулканогенных осадков.

Рис. 4. Общая схема модели с наложенными расчетными изотермами (слева) и изобарами (справа). Условные обозначения: 1 - sediments (S), 2 - continental crust (CÇ), 3 - subcontinental mantle lithosphère (SML), 4 - oceanic crust (OC), 5 - oceanic mantle lithosphère (OML), 6 - asthenosphere (A).

Изолинии: 7 - температуры (°C), 8 - давления (бар).

Принято, что поры заполнены морской водой. Пористость рассчитана по справочной литературе. Зона ОС по химическому составу отвечает комбинации «океанического базальтового» (79.60%), «океанического вулканогенного» (13.92%) и «океанического осадочного» (3.52%) слоев с обводненными порами (2.96%). Зона S комбинирует в себе вулканогенные (72.42%) и морские (14.73%) осадки и насыщена водой до 12.85 масс. %. Состав морской воды нормирован из литературных данных. Зона OML («oceanic mantle lithosphere») представляет

собой переходную к астеносфере зону «базальтового» слоя океанов и свободной воды не содержит. Десятые и сотые доли в концептуальных составах зон появились вследствие пересчета от объемных процентов к массовым.

Табл. 1. Концептуальный состав принять« в модели зон, масс. %.

Слой А БМЬ сс ОМ1. ОС в Водонасы-щенность

Океанический осадочный - - - - 3.52 14.73 20.00

Океанический вулканогенный - - - - 13.92 72.42 25.00

Океанический базальтовый - - - 100 79.60 - 1.50

Континентальный гранитный - - 100 - - - 0.00

Континентальный базальтовый - 100 - - - - 0.00

Астеносфера (пиролит) 100 - - - - - 0.00

Морская вода - - - - 2.96 12.85 100.00

Физический этап моделирования. В расчетах полей температуры использовалась авторская модификация известного метода Гаусса-Зейделя, включающая точки с фиксированными значениями. Модификация заключается в том, что для некоторой совокупности N членов исходной матрицы (точек плоскости с определенными координатами) задаются постоянные значения, и пересчет этих значений в итерациях не происходит. Совокупность N выпадает из итерационных расчетов, но значения ее членов влияют на окружающие. Влияние постепенно ослабевает с удалением от фиксированной точки. Данное развитие метода Гаусса-Зейделя было реализовано в виде авторского программного продукта \ПасН который позволяет рассчитывать и визуализировать двумерные стационарные распределения интенсивных физических величин (в данном случае температуры) при влиянии частных коэффициентов с учетом фиксированных точек и граничных условий.

В нашем случае фиксированными точками являлись подсистемы поверхностного уровня (4°С и 20°С для дна океана и континента соответственно) и подсистемы нижнего уровня астеносферы (1300°С), подсистемы корового «вулканического очага» (400°С) и другие. Значения температуры для подсистем левой и правой границ модели, а также подсистем нижнего уровня зон ОС и ОМЬ рассчитывались рекурсивно. Температуропроводности зон при стандартных условиях сначала рассчитывались по диаграммам (Теплофизические..., 1987), затем, с использованием уравнения Тихомирова, - для частных температур подсистем. Расчетное распределение температурных полей модели показано на рис. 4 (слева).

Литостатическое давление рассчитывалось, исходя из плотностей подсистем, и в нижней части разреза достигло ~ 30 кбар. Полученные изобары показаны на рис. 4 (справа).

Физико-химический этап моделирования. Наборы и содержания независимых компонентов зон и, следовательно, подсистем модели были сформированы в соответствии с методикой, исходя из принятых составов зон. Не ставя себе задачу исследования параметров распределения редких и рассеянных элементов, мы ограничились выбором 13-ти независимых компонентов.

Таблица 2. Независимые компоненты условных зон модели. Под аббревиатурой зоны показана ее плотность (г/смд).

Компонент А р=3.3 8МЬ р =2.7 СС р =2.5 ОМЬ р =3.2 ОС р =3.0 8 р =2.2

$1, моль/кг 6.68 9.18 10.70 8.25 7.90 6.71

А1, моль/кг 0.73 2.82 3.08 3.05 2.88 2.28

Ре, моль/кг 2.12 1.32 0.80 1.42 1.34 1.03

М^ моль/кг 9.10 1.59 0.75 1.96 1.81 1.26

Са, моль/кг 0.44 1.45 0.69 2.01 2.16 2.76

К, моль/кг 0.04 0.28 0.63 0.05 0.07 0.15

Ыа, моль/кг 0.18 0.76 0.89 0.84 0.82 0.77

Н, моль/кг 0.22 1.56 0.78 0.77 4.47 17.00

О, моль/кг 26.59 28.64 29.70 27.84 28.65 31.34

в, моль/кг 0.0000 0.0281 0.0327 0.0181 0.0176 0.0156

С, моль/кг 0.0002 0.0572 0.1313 0.0000 0.1840 0.8759

С1, моль/кг 0.0001 0.0020 0.0060 0.0000 0.0161 0.0699

V, моль/кг 0.0001 0.0089 0.0276 0.0095 0.0106 0.0147

Рис. 5. Расчетное распределение фазы раствора, газовой фазы и свободного флюида в моделируемой метасистеме, масс. %.

Раствор

О 10 2 0 30 40 Расстояние, км

Расстояние, юл

В их состав вошли элементы оксидов стандартного силикатного анализа: 81, А1, Ре, Са, Ыа, К, Н и О, дополненные летучими: 8, С, С1 и Р (табл. 2). Зависимые компоненты системы были отобраны из баз данных ПК «Селектор». Это 158 возможных конденсированных фаз, 98 компонентов водного раствора и 19 компонентов газовой фазы. Расчет на ПК «Селектор» производился методом минимизации изобарно-изотермического потенциала С?(Т, Р) для более 4000 вариантов. Р-Т-условия и вектор Ь каждого варианта соответствовали геометрическому месту подсистем модели.

По всей мегасистеме в целом расчет показал возможность равновесных ненулевых концентраций 49 конденсированных фаз, 61 компонента водного раствора и 6 компонентов газовой фазы. Для всех фаз и компонентов были построены диаграммы концентраций в модельной плоскости. Расчет показал также, что в 2160 резервуарах в равновесии с конденсированными фазами имеется фаза водного раствора и в 786 - газовая фаза. Чаще всего они взаимоисключают друг друга. Обратим внимание на поля существования свободного флюида в области дегидратации погружающейся плиты на глубинах 50-100 км (рис. 5).

«Сухой» солидус базальта

Флюид, масс.%

4 50 00 3.50 300 2.50 2.00 1 50 1 00

Щц

Йо.50

600 800 1000 1200 Температура. С

Рис. 6. Расчетная Р-Т диаграмма распределения свободного флюида с линиями солидуса и ликвидуса базальта (Перчук, 1973). Квадратами обозначены три условно выделенные зоны повышенного содержания свободного флюида.

На Р-Т-диаграмме можно выделить три области повышенного содержания свободного флюида, которые попадают в два участка частичного и один участок практически полного плавления «мокрого» базальта (рис. 6). Это согласуется с литературными данными (Перчук, 1973; Добрецов и др., 2001).

Динамический этап моделирования. Расчет полной статической картины распределения вещества модели позволил перейти к исследованию физико-химического взаимодействия подсистем между собой для моделирования подъема флюида (группы подвижных фаз).

Для этого автором использовался метод последовательного проточного реактора с 20 резервуарами, расположенными снизу вверх в вертикальном направлении вдоль правой границы метасистемы.

Группа подвижных фаз (флюид) включала газовую фазу и раствор. Геометрическая область максимальных концентраций свободного флюида в условиях дегидратации пород океанической плиты соответствует на Р-Т диаграмме ~800°С и 30360 бар. Эти условия приняты для 1 -го резервуара.

За состав исходного флюида принят расчетный равновесный свободный флюид, соответствующий условиям погружения плиты на глубину 95-100 км (резервуар 1). Параметры последнего 20-го резервуара (глубина 5 км) составили 84°С и 1250 бар. Характеристики принятых в модели подъема флюида резервуаров приведены в табл. 3.

Таблица 3. Параметры модельных резервуаров при расчете подъема флюида.

№ резервуара Глубина, км Т,°С Р, бар Зона

20 5 84 1250 Continental crust (СС)

19 11 256 2750

18 16 376 4000

17 21 441 5250

16 26 508 6500

15 31 578 7770 Subcontinental mantle lithosphere (SML)

14 37 661 9390

13 43 744 11070 Asthenosphere (A)

12 48 814 12720

11 53 886 14370

10 58 959 16020

9 63 1024 17670

8 68 1077 19320

7 73 1131 20970

6 78 1115 22620

5 83 1052 24270

4 88 986 25920

3 93 908 27570

2 97 846 28860 Oceanic crust (ОС)

1 100 786 30360

На рис. 7 приведен расчетный характер изменения содержаний компонентов флюида на разных глубинах относительно от его первоначального состава на глубине 100 км.

Рис. 7. Изменения концентраций независимых компонентов флюида, нормированные к исходному составу. Шкала концентраций - логарифмическая.

Из диаграммы очевидна тенденция к увеличению концентраций большинства компонентов по мере подъема флюида в верхние горизонты, что наиболее характерно для натрия, магния и железа, в меньшей степени - для кальция, алюминия, серы и галогенов.

Исключение составляют кремний, наиболее интенсивно расходующийся на поддержание силикатных фаз уже в астеносферном клине и в меньших масштабах - калий, но уже на уровне коровых резервуаров. Также важным результатом моделирования эволюции флюида явилось расчетное отсутствие равновесных фаз раствора и газа в астеносферном клине в интервале глубин 60-85 км, где компоненты флюида конденсируются в твердые минеральные фазы или расплав.

Расчетные минеральные фазы по характеру изменений и приуроченности к конкретным резервуарам можно весьма определенно разделить на четыре группы. Критерием такого разделения выступают, во-первых, геометрическое место (интервалы) нахождения измененных и неизмененных минеральных фаз, а во-вторых, количественные показатели изменения содержания того или иного минерала. В первую группу входят минералы, вариации содержаний которых приурочены к верхним частям литосферы. В качестве примера минералов этой группы приведем кварц, калиевый полевой шпат, гематит, энстатит. Амплитуда изменений содержаний минералов этой группы относительно высока и лежит в пределах 8-16 масс. %. Вторая группа минеральных фаз, напро-

тив, претерпевает значительные (до 8-16 масс. %) изменения в «нижних» резервуарах (до 14-го включительно). Их содержание падает в 3-4 резервуарах, затем колеблется до 14-го резервуара, после которого изменения прекращаются. Это такие минеральные фазы как периклаз, оливин, вюстит и магнетит. Третья и четвертая группы, в отличие от первых двух, характеризуются невысокими амплитудами вариаций составов (менее 1 масс. %), но качественно повторяют их поведение. Минералы третьей группы изменяются в резервуарах выше 14-го, четвертой группы - в нижележащих резервуарах.

Наибольший интерес вызывает изменение в поведении практически всех конденсированных фаз на рубеже 14-го резервуара, который в модели соответствует глубине 35 км и переходу «астеносфера -нижняя кора». Причиной этому, несомненно, служит резкое различие в составах модельных зон «А» и «SML», а также значительный градиент температуры на данном участке вертикального движения флюида.

Основные выводы по модели зоны субдукиии. Для таких длительно живущих объектов, как зоны субдукции, распределение интенсивных параметров можно с достаточной степенью уверенности считать стационарным. Поэтому применение расчетного стационарного распределения температурных полей на основании совокупности точек с эмпирически или теоретически известными значениями, на наш взгляд, вполне приемлемо. Рассчитанные значения качественно соответствуют общепринятой температурной модели (Honda, Uyeda, 1983; Добрецов и др., 2001) с несколько повышенными значениями для субдуцирующей плиты и пониженными - для астеносферного клина.

Основанные на полученном распределении температур физико-химические расчеты подтверждают классическое разбиение погружающейся плиты на сегменты (Добрецов, Кирдяшкин, 1997), несколько сместив их границы. Верхний (<30 км) характеризуется условиями фации зеленых сланцев при температуре меньше 400°С с незначительной дегидратацией. Следующий сегмент (30-50 км) представляет собой область фации глаукофановых сланцев при 400-700°С. Третий сегмент (>50 км) характеризуется температурами 700-1000°С эклогитовой фации с интенсивной дегидратацией и частичным плавлением. Самая нижняя (~100 км) часть рассматриваемой зоны захватывает самое начало четвертого сегмента - области более полного плавления плиты — возможного источника расплава для вулканизма II типа по классификации (Добрецов и др., 2001).

Расчеты также показали невозможность прохождения флюида через астеносферный клин, в резервуарах которого нет подвижных фаз. Растворы коры, таким образом, имеют астеносферный генезис.

Заключение

После многолетней (1997-2009) работы в данной области, в течение которой были разработаны, реализованы и апробированы различные прикладные компьютерные модели, автору удалось интегрировать полученный опыт в единую методику комплексного компьютерного моделирования. За эти годы были разработаны разнотипные модели различных геохимических объектов, написаны и протестированы необходимые прикладные программные продукты: численная конечно-разностная модель тепло- и массопереноса в гидротермальной системе СОХ; численная термодинамическая модель захоронения токсичных отходов в паровых зонах гидротермальных систем областей развития современного вулканизма; численная модель формирования гидротермального золоторудного месторождения; программное обеспечение для расчета AFC-моделей эволюции расплава; статистическая компьютерная модель эволюции флюидозаполненных трещин в пластичной среде; компьютерная геодинамическая модель коллизии плит; численная физико-химическая модель образования рудоносных гранитоидов; комплексная компьютерная модель зоны субдукции и др.

Цель работы - разработка методики комплексного компьютерного моделирования геохимических объектов на основе объектно-ориентированного подхода и апробация ее на прикладных тематических моделях геологических объектов - на данный момент выполнена в полном объеме, как и поставленные задачи. Автором весьма детально освоена теория компьютерного моделирования, а также языки программирования Basic, Pascal, С, С++ и соответствующие интегрированные среды разработки, что позволило создавать собственные модели и программные продукты. Были детально проанализированы существующие методы физико-химического моделирования. Во время работы автор принимал активное участие в тестировании и апробации различных версий программного комплекса «Селектор» как в Лаборатории физико-химического моделирования Института геохимии СО РАН (г. Иркутск), так и в Лаборатории геохимии Геологического института СО РАН (г. Улан-Удэ).

Результаты применения предлагаемой методики в комплексном компьютерном моделировании зон субдукции, магматогенно-гидротермальных систем и коллизии плит позволяют подтвердить некоторые теоретические догадки для сложных в изучении глубоких горизонтов данных объектов и поднять будущие модельные работы на новый, детальный уровень. Подведена вещественная основа для дальнейших исследований в области зон субдукции в виде базы данных, включающей более четырех тысяч равновесных минеральных ассоциаций с

раствором и газовой фазой, рассчитанных исходя из наиболее общих исходных данных. Также заслуживает внимание расчетное распределение температурных полей, наиболее приемлемое с точки зрения петро-физики геологических сред.

Непосредственное развитие своих исследований автор видит в создании нового программного комплекса для расчетов гидродинамической эволюции растворов в трещиновато-пористых средах, органичном соединении этого программного комплекса с ПК «Селектор» и в создании и расчете новых комплексных компьютерных моделей различных геохимических объектов с постоянным повышением вышеописанных параметров корректности.

СПИСОК ОСНОВНЫХ ОПУБЛИКОВАННЫХ РАБОТ ПО ТЕМЕ ДИССЕРТАЦИИ

1. Васильев В.И., Жатнуев Н.С. Физико-химическое моделирование минеральных равновесий зон субдукции островных дуг на программном комплексе «Селектор». // Вулканизм и геодинамика: Материалы III Всероссийского симпозиума по вулканологии и палеовулканологии. - Т.1. - Улан-Удэ: Изд-во БНЦ СО РАН, 2006. - С. 128132.

2. Васильев В.И. Численное моделирование минеральных равновесий зон субдукции на ПК «Селектор». // Геохимия и рудообразование радиоактивных, благородных и редких металлов в эндогенных и экзогенных процессах. Материалы Всероссийской конференции с иностранным участием. Часть 2. - Улан-Удэ: Изд-во БНЦ СО РАН, 2007.-С. 116-118.

3. Васильев В.И., Жатнуев Н.С. Реализация модели распределения вещества и тепла при коллизии на языке СИ++ с привлечением ПК «Селектор». // Геохимия и рудообразование радиоактивных, благородных и редких металлов в эндогенных и экзогенных процессах. Материалы Всероссийской конференции с иностранным участием. Часть 2. - Улан-Удэ: Изд-во БНЦ СО РАН, 2007. - С. 119-121.

4. Васильева Е. В., Васильев В.И., Жатнуев Н.С. Тектонофизическое моделирование динамики флюидсодержащих трещин в литосфере. // Геохимия и рудообразование радиоактивных, благородных и редких металлов в эндогенных и экзогенных процессах. Материалы Всероссийской конференции с иностранным участием. Часть 2. -Улан-Уда: Изд-во БНЦ СО РАН, 2007. -С. 122-125.

5. Vasiliev V.l., Khrustalev V.K. The numerical thermodynamic model of ore-bearing PZ2.r-granitoids of Vitim plateau. II Granites and Earth's Evolution: Geodynamic setting, Pedogenesis and Ore Content of Granitoid Batholiths: Proceedings of the First International Geologic Conference. - Ulan-Ude: Publishing House BSC SB RAS, 2008. P. 54-56.

6. Васильев В.И. Численная физико-химическая модель подъема свободного флюида в зоне субдукции. // Современные проблемы геологии, геохимии и геоэкологии Дальнего Востока России. Материалы II региональной конференции- Владивосток: Дальнаука, 2008.-С.108-110.

7. Васильев В.И., Чудненко К.В., Жатнуев Н.С., Васильева Е.В. Комплексное компьютерное моделирование геологических объектов на примере разреза зоны субдукции. // Геоинформатика, №3,2009. - С. 15-30.

8. Васильев В.И. Массоперенос и минералообразование в магматогенно-гидротермальных системах по результатам численного физико-химического моделирования. // Вулканизм и геодинамика: Материалы IV Всероссийского симпозиума по вулканологии и палеовулканологии. - Т. 2. - Петропавловск-Камчатский: ИВиС ДВО РАН, 2009. - С. 713-717.

Подписано в печать 16.10.2009 г. формат 60x84 1/16. Бумага офсетная. Усл. печ. л.1,3 Тираж 100. Заказ № 47

Отпечатано в типографии Изд-ва БНЦ СО РАН. 670047 г. Улан-Удэ ул. Сахъяновой, 6.

Содержание диссертации, кандидата геолого-минералогических наук, Васильев, Владимир Игоревич

Введение

Глава 1. Состояние проблемы

1.1. Геометрический аспект моделирования

1.2. Физический аспект моделирования

1.3. Физико-химический (термодинамический) аспект моделирования

1.4. Некоторые классические модели

1.4.1. Модель охлаждения интрузива JI.M. Кэтлса

1.4.2. Модели тепломассопереноса A.B. Кирюхина и В.М. Сугробова

1.4.3. Комбинированные модели Д.В. Гричука,

М.В. Борисова и A.B. Тутубалина

Глава 2. Методика моделирования

2.1. Основные определения

2.2. Геометрический этап

2.3. Физический этап

2.4. Физико-химический этап

2.5. Динамический этап

Глава 3. Модель-схема коллизии плит

3.1. Описание модели —

3.2. Программный продукт Vladi Collision 2.

Глава 4. Модель магматогенно-гидротермальной системы

Глава 5. Модель зоны субдукции

5.1. Концептуальная модель —

5.1.1. Постановка задачи —

5.1.2. Петрогеохимическая характеристика островных дуг Курило-Камчатской системы

5.1.3. Природа проявления современного вулканизма Курило-Камчатской системы

5.1.4. Модель магмообразования в зоне субдукции

5.2. Геометрический анализ и входные составы зон модели

5.3. Физический этап

5.3.1. Метод Гаусса-Зейделя с фиксированными точками —

5.3.2. Программный продукт У1асН 1.

5.3.3. Расчет температур и давлений подсистем модели

5.4. Физико-химический этап

5.5. Динамический этап

Введение Диссертация по наукам о земле, на тему "Комплексное компьютерное моделирование геохимических объектов на примерах двумерных моделей коллизии плит, магматогенно-гидротермальной системы и зоны субдукции"

Актуальность проблемы. Общей тенденцией развития современной геохимии является переход от преимущественно описательной характеристики геологических объектов к изучению тех процессов, которые управляют превращениями вещества в недрах и на поверхности нашей планеты. Научную ценность химико-аналитического материала, количественным показателям которого еще недавно придавали первостепенное значение, предпочитают рассматривать с позиций того, что нового может дать этот материал в раскрытии физических и физико-химических закономерностей природных процессов.

Смещение фундаментальных научных интересов от сбора и накопления аналитической информации к исследованию процессов преобразования, миграции и эволюции вещества на теоретическом уровне строгих положений химической термодинамики объективно продиктовано естественным ходом развития геохимии и геофизики как наиболее точных областей человеческого знания в комплексе наук о Земле. В то же время использованием компьютера в качестве важнейшего инструмента исследований в этой области пренебрегает большинство ученых.

Между тем, дальнейший рост научного потенциала геохимической науки зависит уже не столько от тех дорогостоящих усилий, которые направлены на аддитивное увеличение общей суммы аналитических данных, сколько от глубины и строгости физико-химического синтеза всей многосложной совокупности полевой, аналитической и экспериментальной информации в науках о Земле и смежных областях знания. Такого уровня исследований на современном этапе можно добиться только методами компьютерного моделирования.

Подводя итог вышесказанному, подчеркнем насущную необходимость универсального методологического подхода к компьютерному моделированию природных объектов, способного стать инструментом изучения уровня любого масштаба геохимических исследований. Также исключительно важна популяризация компьютерного моделирования среди геологов различной специализации, демонстрация его возможностей на примерах тематически близких им моделей.

Целью работы является разработка методики комплексного компьютерного моделирования геологических объектов на основе объектно-ориентированного подхода и апробация ее на прикладных тематических моделях геологических объектов. Это модель-схема коллизии плит, физико-химическая модель магматогенно-гидротермальной системы и многорезерву-арная модель обстановки зоны субдукции.

Задачи работы были поставлены в соответствии с выбранной целью и могут быть сформулированы следующим образом:

1) освоение теории компьютерного моделирования;

2) освоение языков программирования и интегрированных сред разработки (языки Basic, Pascal, С, С++);

3) освоение и модернизация существующих программных продуктов для различных типов моделирования геологических объектов; анализ существующих методов физико-химического моделирования;

4) наработка опыта в построении существенного количества разнотипных компьютерных моделей; разработка авторских прикладных программ, необходимых при построении моделей разного типа;

5) собственно разработка и апробация методики;

6) анализ, интерпретация и оценка результатов моделирования.

Научная новизна. Существует множество методов компьютерного моделирования с частных позиций геометрического, физического, химического или физико-химического подходов, иногда парно объединенных. Тем не менее, попытка свести разноплановые методы в схему единого универсального подхода к моделированию геологических объектов, выполнена, по-видимому, впервые. Кроме этого, по мнению автора, опубликованные в доступной научной литературе модели являются результатом такого положения вещей, когда есть множество методов, но нет методик компьютерного моделирования, т.е. отсутствует единый методологический принцип разработки и оценки моделей геологических объектов. В данной работе автор постарался восполнить этот недостаток современного этапа исследования.

Практика применения комплексного подхода в компьютерном моделировании, включающего аспекты представления и описания исследуемого объекта в рамках перехода от эмпирической модели к компьютерной, получила достаточное развитие в первую очередь в приложениях к решению инженерно-конструкторских задач (Нигматулин, 1987; Оран, Борис, 1990; Седов, 2004). Высокая степень детерминированности изучаемых объектов и детально разработанный математический аппарат делает такие реализации своего рода эталонами прикладных методов математического моделирования. В моделировании же природных систем далеко не всегда могут быть использованы передовые разработки из области технологических приложений, так как геохимические процессы формирования и развития реальных геологических объектов обусловлены протеканием многостадийных химических превращений вещества, кинетические характеристики которых не всегда известны с достаточной степенью надежности и достоверности; поэтому построение детерминированной модели в рамках изначально заданной системы дифференциальных уравнений с учетом всего многообразия обратимых и необратимых химических взаимодействий на существующем уровне познания не представляется возможным.

В то же время наибольшие успехи в построении физико-химических моделей природных процессов получены с помощью методов термодинамического моделирования (Карпов, 1981), которые позволяют производить расчет компонентного и фазового состава моделируемых систем без рассмотрения траекторий перехода систем в конечное равновесное состояние, опираясь на фундаментальный принцип достижения системы минимума термодинамического потенциала. Дальнейшее развитие методов термодинамического моделирования в геологии направлено на преодоление главного ограничение аппарата химической термодинамики — отсутствие пространственно-временных координат и связано с разработкой подходов, объединяющих динамические и физико-химические модели в рамках единых вычислительных схем (Шарапов и др., 2007, 2008).

Объектно-ориентированный подход, широко применяемый в программировании, при моделировании природных объектов также является новым; работы такого плана автору не известны. Между тем, именно он позволяет максимально наглядно и удобно для восприятия представлять модели любых геологических объектов в терминах конструктивных схем единой концепции.

Нельзя не отметить также несомненную новизну компьютерных моделей (главы 3-5), разработанных на базе оригинальной методики, и программных продуктов, разработанных автором по собственным алгоритмам.

Практическая значимость работы состоит в универсальности применения представляемой методики для самого широкого круга геологических задач, в которых в той или иной степени затрагивается проблема эволюции вещества. Это утверждение аргументируется представлением методологически единых моделей таких категориально различных геологических объектов, как процесс коллизии плит и явление термодинамического резервуара в толще породы (глава 2). Разработанные модели (главы 3-5) показывают практическую пригодность методики для получения как качественных, так и количественных параметров моделируемого объекта.

Результаты применения предлагаемой методики в области комплексного моделирования зоны субдукции позволяют не только подтвердить многие теоретические догадки для сложных в изучении глубоких горизонтов данных областей, но и поднять будущие модельные работы на новый, детальный уровень. Так, подведена серьезная вещественная основа для дальнейших исследований в этой области: создана база данных, включающая более четырех тысяч равновесных минеральных ассоциаций с раствором и газовой фазой, рассчитанных исходя из наиболее общих исходных данных. Также заслуживает внимание расчетное распределение температурных полей, наиболее приемлемое с точки зрения петрофизики геологических сред.

Кроме собственно геолого-геохимической практической значимости, заслуживают внимания предлагаемые относительные критерии оценки корректности моделей (глава 2). Автор убежден, что будущие модельные исследования необходимо вести в направлении роста именно этих критериев, оценивая любую компьютерную модель хотя бы сравнительно с другими моделями относительно степени учета основных характеристик и функциональности на примере подобных геологических объектов.

Защищаемые положения.

1. Предложенная методика комплексного компьютерного моделирования позволяет в рамках универсального объектно-ориентированного подхода разрабатывать модели самых различных геохимических объектов и оценивать их по предложенным критериям корректности: описательному, функциональному и событийному.

2. Модель коллизии плит, разработанная по предложенной методике и реализованная с помощью прикладного программного продукта Vladi Collision, позволяет, варьируя начальными и граничными условиями, рассчитывать и визуализировать различные варианты геометрической эволюции модельной системы с учетом распределения температурных полей разреза.

3. Созданная модель магматогенно-гидротермальной системы отражает основные закономерности массопереноса и минералообразования в гидротермальных системах областей развития современного вулканизма.

4. Разработанная комплексная компьютерная модель зоны субдукции, представленной в виде двумерной метасистемы, за счет разбиения на -4000 подсистем, подтверждает гипотетические положения о Р-Т—условиях нахождения свободного флюида при его подъеме от дегидратирующейся плиты к поверхности.

Реализация работы. Различные варианты авторских программных продуктов переданы и используются в лабораториях Геологического института СО РАН (г. Улан-Удэ) и Института геохимии СО РАН (г. Иркутск).

Публикации и апробация. Основные положения диссертации опубликованы в восемнадцати печатных работах, а также докладывались на Школе-семинаре российских делегатов XXXI Международного Геологического Конгресса (БИС «Академик Иоффе», Атлантика, 2000), XIX Всероссийской конференции «Строение литосферы и геодинамика» (Иркутск, 2001), III Всероссийском симпозиуме с международным участием «Золото Сибири и Дальнего Востока» (Улан-Удэ, 2004), III Всероссийском симпозиуме по вулканологии и палеовулканологии (Улан-Удэ, 2006), Всероссийской конференции с иностранным участием в честь 50-летия СО РАН и 80-летия чл.-корр. РАН Ф.П. Кренделева (Улан-Удэ, 2007), I международной конференции «Граниты и эволюция Земли» (Улан-Удэ, 2008), IV Всероссийском симпозиуме по вулканологии и палеовулканологии «Вулканизм и геодинамика» (Петропавловск-Камчатский, 2009), Научных совещаниях «Геодинамическая эволюция литосферы Центрально-Азиатского подвижного пояса (от океана к континенту)» (Иркутск, 2008, 2009), ежегодных Научных сессиях Геологического института СО РАН (1998 - 2009).

Объем и структура работы. Диссертация состоит из введения, пяти глав, заключения и списка литературы общим объемом 120 страниц печатного текста, 16 таблиц и 36 иллюстраций.

Заключение Диссертация по теме "Геохимия, геохимические методы поисков полезных ископаемых", Васильев, Владимир Игоревич

5.6. Основные выводы

В представленной модели были рассмотрены все четыре основных аспекта моделирования геологического объекта: геометрический, физический, физико-химический и динамический.

За счет разбиения метасистемы на большое количество подсистем и учета доступного количества их экстенсивных параметров, на геометрическом этапе была достигнута доскональная проработанность структуры модели. Наименьшая по размерам зона (5) включала в себя 48 подсистем, наибольшая - 1453 (ОМЬ). Такой подробный геометрический анализ позволил детально рассчитать распределение тепла и вещества на физическом и физико-химическом этапах.

Для таких длительно живущих геологических объектов, как зоны субдук-ции, распределение интенсивных параметров можно с достаточной степенью уверенности считать стационарным. Поэтому применение расчетного стационарного распределения температурных полей на основании совокупности точек с эмпирически или теоретически известными значениями, на наш взгляд, вполне приемлемо. В целом, рассчитанные значения температур качественно соответствуют общепринятой температурной модели (Honda, Uyeda, 1983; Добре-цов и др., 2001) с несколько повышенными значениями для субдуцирующей плиты и пониженными - для астеносферного клина.

Основанные на полученном распределении температур физико-химические расчеты позволяют подтвердить классическое разбиение зоны суб-дукции плиты на сегменты (Добрецов, Кирдяшкин, 1997), несколько сместив их границы. Верхний сегмент (<30 км) характеризуется условиями фации зеленых сланцев при температуре меньше 400°С с незначительной дегидратацией. Следующий сегмент (30—50 км) представляет собой область фации глаукофановых сланцев при 400-700°С. Третий сегмент (>50 км) характеризуется температурами 700-1000°С эклогитовой фации с интенсивной дегидратацией и частичным плавлением. Самая нижняя (-100 км) часть рассматриваемой зоны захватывает самое начало четвертого сегмента - области более полного плавления погружающейся плиты — возможного источника расплава для вулканизма II типа по классификации (Tatsumi, 1989; Добрецов и др., 2001).

Результаты моделирования позволяют не только подтвердить многие теоретические догадки для сложных в изучении горизонтов зон субдукции, но и поднять будущие модельные работы на новый, детальный уровень. Так, подведена серьезная вещественная основа для дальнейших исследований в этой области: создана база данных, включающая более четырех тысяч равновесных минеральных ассоциаций с раствором и газовой фазой, рассчитанных исходя из наиболее общих исходных данных. Также заслуживает внимание расчетное распределение температурных полей, наиболее обоснованное математически и теплофизически (Теплофизические., 1987; Васильев, Жатнуев, 2006; Васильев и др., 2009).

В результате моделирования было рассчитано изменение основных характеристик флюида при его подъеме от дегидратирующейся плиты к поверхности через цепочку последовательно сопряженных транзитных подсистем (резервуаров). Расчеты показали невозможность прямого прохождения флюида через ас-теносферный клин, в резервуарах которого подвижные фазы отсутствуют. Растворы нижней коры, таким образом, имеют явно выраженный астеносферный генезис.

ЗАКЛЮЧЕНИЕ

После многолетней (1997-2009) работы в данной области, в продолжение которой были разработаны, реализованы и апробированы различные прикладные компьютерные модели, автору удалось интегрировать полученный опыт в единую методику комплексного компьютерного моделирования. За эти годы были разработаны разнотипные модели различных геологических объектов, написаны и протестированы необходимые прикладные программные продукты. Это численная конечно-разностная модель тепло- и массопереноса в гидротермальной системе СОХ (Васильев, 2000, 2001); численная термодинамическая модель захоронения токсичных отходов в паровых зонах гидротермальных систем областей развития современного вулканизма (Отчет по проекту №401 VI Конкурса-экспертизы молодых ученых РАН по фундаментальным и прикладным исследованиям, 2001); численная модель формирования гидротермального золоторудного месторождения (Калинин и др., 2004); физико-химическая модель формирования золотоносных кор выветривания Восточного Саяна (Калинин и др., 2005; Жмодик и др., 2005); модель эволюции магматического очага и подъема магмы по цилиндрическому каналу (Жатнуев и др., 2006); программное обеспечение для расчета АРС-моделей эволюции расплава в магматической камере (Отчет по проекту №120 Лаврентьевского конкурса СО РАН, 2006); статистическая компьютерная модель эволюции флюидозаполненных трещин в пластичной среде (Васильева, Васильев, Жатнуев, 2007); компьютерная геодинамическая модель коллизии плит (Васильев, Жатнуев, 2007); численная физико-химическая модель образования рудоносных гранитоидов (УавШеу, Юиг^а-1еу, 2008); комплексная компьютерная модель зоны субдукции (Васильев, 2006, 2007) и другие.

Цель диссертационной работы — разработка методики комплексного компьютерного моделирования геологических объектов на основе объектно-ориентированного подхода и апробация ее на прикладных тематических моделях геологических объектов - на данный момент выполнена в полном объеме, как и поставленные задачи. Автором весьма детально освоена теория компьютерного моделирования, а также языки программирования Basic, Pascal, С, С++ и соответствующие интегрированные среды разработки, что позволило не только создавать собственные модели и программные продукты, но и на протяжении восьми лет (2000—2008) преподавать на кафедре геологии Бурятского государственного университета такие предметы как «Информатика», «Компьютерное моделирование», «ЭВТ в геологии» и «Основы физико-химического моделирования».

Были детально проанализированы существующие методы физико-химического моделирования. Изначальная приверженность к вышеописанному «методу реакций» постепенно сменилась твердой убежденностью в правильности «метода минимизации». Во время работы автор принимал активное участие в тестировании и апробации различных версий программного комплекса «Селектор» как в Лаборатории физико-химического моделирования Института геохимии СО РАН (г. Иркутск), так и в Лаборатории геохимии Геологического института СО РАН (г. Улан-Удэ).

Представляемая методика основана на универсальном объектно-ориентированном подходе и, несомненно, заслуживает не только права на существование и развитие, но и популяризации среди исследователей в области наук о Земле.

Результаты применения предлагаемой методики в комплексном компьютерном моделировании зон субдукции, магматогенно-гидротермальных систем и коллизии плит позволяют не только подтвердить многие теоретические догадки для сложных в изучении глубоких горизонтов данных объектов, но и поднять будущие модельные работы на новый, детальный уровень. Так, подведена серьезная вещественная основа для дальнейших исследований в области зон субдукции: создана база данных, включающая более четырех тысяч равновесных минеральных ассоциаций с раствором и газовой фазой, рассчитанных исходя из наиболее общих исходных данных. Также заслуживает внимание расчетное распределение температурных полей, наиболее приемлемое с точки зрения петрофизики геологических сред.

Непосредственное развитие своих исследований автор видит в создании нового программного комплекса для расчетов гидродинамической эволюции растворов в трещиновато-пористых средах, органичном соединении этого программного комплекса с ПК «Селектор» и в создании и расчете новых комплексных компьютерных моделей различных геологических объектов с постоянным повышением вышеописанных параметров корректности.

Библиография Диссертация по наукам о земле, кандидата геолого-минералогических наук, Васильев, Владимир Игоревич, Улан-Удэ

1. Авдейко Г.П. Геодинамика проявления вулканизма Курильской островной дуги и оценка моделей магмообразования. // Геотектоника, №2, 1994. -С. 19-32.

2. Авдейко Г.П., Волынец О.Н., Антонов А.Ю. Вулканизм Курильской островной дуги: структурно-петрологические аспекты и проблемы магмообразования. // Вулканология и сейсмология, №5, 1989. С. 3-16.

3. Авдейко Г.П., Пилипенко Г.Ф., Палуева A.A., Напылова O.A. Геотектонические позиции современных гидротермальных проявлений Камчатки. // Вулканология и сейсмология, №6, 1998. С. 85-99.

4. Авдейко Г.П., Попруженко C.B., Палуева A.A. Современная тектоническая структура Курило-Камчатского региона и условия магмообразования. / Геодинамика и вулканизм Курило-Камчатской островодужной системы. -Петропавловск-Камчатский, 2001. С. 9-34.

5. Авдейко Г.П., Попруженко C.B., Палуева A.A. Тектоническое развитие и вулкано-тектоническое районирование Курило-Камчатской островодужной системы. //Геотектоника, №4, 2002. — С. 64-80.

6. Аверкин И.А., Шарапов В.Н. Динамика массообмена в гидротермальных системах на температурном геохимическом барьере. Новосибирск: Изд-во ИГИГ, 1986.-46 с.

7. Аплонов C.B. Геодинамика. СПб: Изд-во СПбГУ, 2001. - 360 с.

8. Богатиков O.A., Косарева JI.B., Шарков Е.В. Средние химические составы магматических горных пород. -М.: Недра, 1987. 152 с.

9. Большое трещинное Толбачинское извержение (1975 — 1976 гг., Камчатка). -М.: Наука, 1984. 638 с.

10. Борисов М.В., Шваров Ю.В. Термодинамика геохимических процессов. — М.: Изд-во МГУ, 1992. 256 с.

11. Быков А. Желаемое и действительное в геометрическом моделировании. // САПР и графика, №1, 2002. С. 3 - 26.

12. Васильев В.И. Отчет по проекту №401 VI Конкурса-экспертизы молодых ученых РАН по фундаментальным и прикладным исследованиям, 2001.

13. Васильев В.И. Численное моделирование динамики тепломассопотоков и минералообразования в гидротермальной системе. // Материалы XIX

14. Всероссийской конференции «Строение литосферы и геодинамика». — Иркутск, 2000. С. 9.

15. Васильев В.К, Хубанов В.Б. Отчет по проекту №120 Лаврентьевского конкурса СО РАН, 2006.

16. Васильев В.И., Чудненко К.В., Жатнуев Н.С., Васильева Е.В. Комплексное компьютерное моделирование геологических объектов на примере разреза зоны субдукции. // Геоинформатика, №3, 2009. С. 15-30.

17. Виноградов А.П. Введение в геохимию океана. М.: Наука, 1967. — 214 с.23.