ТЕРМОДИНАМИЧЕСКИЕ МОДЕЛИ СУБМАРИННЫХ ГИДРОТЕРМАЛЬНЫХ СИСТЕМ

2.2. Моделирование метасоматических процессов в протяженных гидротермальных системах.

    Протяженные гидротермальные системы с трещинной фильтрацией по характеру протекающих в них процессов сильно отличаются от рассматривавшегося выше процесса формирования инфильтрационной метасоматической колонки. Эти отличия выражаются прежде всего в том, что для объемных систем не применимы условия I и II, использованные Д.С.Коржинским. Действительно, при трещинной фильтрации в существенно неоднородных по проницаемости породах возникает значительная гидродинамическая дисперсия, и поршневое вытеснение (условие I) не реализуется. Скорость движения растворов по трещинам может существенно превышать скорости их проникновения внутрь блоков породы и метасоматического изменения внутренних частей блоков. Результатом этого является присутствие останцов неизмененных пород и концентрической зональности метасоматических изменений в блоках породы (нарушение условия II). Вследствие этого метасоматические процессы в таких гидротермальных системах являются по сути комбинацией инфильтрационного метасоматоза в направлении движения раствора и диффузионного метасоматоза по нормали к нему, осложненной сильной гидродинамической дисперсией. Они, очевидно, не могут быть сведены к формированию идеальной инфильтрационно-метасоматической колонки. Строгий анализ таких процессов представляет значительные трудности, и при необходимости решения задач моделирования подобных геологических процессов приходится идти по пути аппроксимаций.

.рис.2.7

    Характерным типом задач моделирования процессов в протяженных системах является модель формирования рудоносных растворов в конвективных гидротермальных системах с экзогенным источником раствора (атмосферным или морским) и эндогенным источником энергии. В этом случае рудоносность раствора определяется взаимодействием вода-порода в термически неоднородной среде с преобладающей трещинной фильтрацией. Наиболее изученный пример таких задач - это модели гидротермальных систем в океанической коре (см. раздел 3.4.1), аналогичный характер имеют модели эпитермальных и стратиформных месторождений [Heinrich et al, 1995; Ilchik, Barton, 1997], и даже грейзенового процесса [Коротаев и др., 1992].

    Для моделирования таких систем нами был предложен метод многоволнового проточного ступенчатого реактора (МПСР) [Гричук, Борисов, 1983; Гричук, 1988]. Конвективная гидротермальная система представляется в нем как политермический ступенчатый реактор (см.рис.2.7, аналог - рис.2.1г), на каждой ступени которого достигается частичное равновесие раствора с породой. Путь течения раствора разбивается на ступени с фиксированным шагом по температуре (10-50^o). На каждой ступени должен быть задан параметр частичного равновесия, определяющий необратимое взаимодействие раствора с породой - "отношение порода/вода" (П/В). Фактически - это критерий подобия модели, представляющий собой отношение массы свежей породы, прореагировавшей с данной порцией раствора на данной ступени к массе этой порции.

Масса порции воды определяется выражением:

, (2.7)

где  - эффективное сечение, v - скорость движения раствора, _w - плотность раствора, t - время прохождения через ступень реактора - температурную зону на пути конвекции, Q - расход раствора. Масса прореагировавшей породы определяется выражением:

,(2.8)

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

    Параметры, входящие в выражения (2.7) и (2.8) возможно оценить сейчас только с очень большими интервалами неопределенности, в особенности - v_r , и они не могут еще служить основой для построения количественных моделей [Рафальский, 1993]. Поэтому в [Гричук, Борисов, 1983; Гричук, 1988] был предложен иной, обходной путь. Дело в том, что для определения величины П/В есть косвенные геохимические методы, в том числе и изотопные [Spooner et al., 1977; Norton, Knight, 1977; и др.]. В частности, для современных океанских гидротермальных систем разработан метод оценки суммарного отношения П/В, достигнутого в гидротермальной системе в целом, по концентрациям легкоподвижных элементов (подробнее см. раздел 3.1.2). Для изученных гидротермальных систем эти величины составляют от 0,5 до 2,0 кг породы на кг раствора. Недостатком таких оценок с точки зрения метода МПСР является то, что величину П/В требуется разделить в модели на вклады каждой ступени. Кроме того, П/В для системы может меняться во времени.

    Для модели конвективной гидротермальной системы можно сделать следующий ряд упрощающих допущений: а) модель стационарна в отношении положения температурных зон и скоростей перемещения раствора; б) времена прохождения температурных зон примерно одинаковы, и время прохождения раствором i-й ступени t_i пропорционально T_i; в) фильтрационные свойства пород меняются во времени и пространстве несущественно. В рамках этих допущений параметры Q и S_fr могут быть приняты одинаковыми по ступеням реактора и не зависящими от времени (хотя и не известными численно). Единственным параметром, входящим в выражение для П/В и зависящим от температуры и времени остается в этом случае удельная скорость реакции со свежей породой v_r .

    Зависимость П/В от температуры в рамках принятых приближений можно записать как

, (2.9)

рис.2.8

где T_o и T_f - температуры начала и конца процесса в нисходящей ветви, T - температурный шаг ступеней реактора, v_ri - средняя скорость реакции в интервале температур, соответствующем ступени i, n - число ступеней в реакторе. Для описания зависимости скорости реакции раствора с породой от температуры можно использовать эмпирическое уравнение, полученное Б.Вудом и Дж.Уолтером [Wood, Walther, 1983; Уолтер, Вуд, 1989], которые показали, что для условий, далеких от равновесия в околонейтральных растворах скорость растворения алюмосиликатов, нормированная на число кислородных атомов в их формулах подчиняется уравнению аррениусовского типа:

, (г-ат О/см^{2 _.}с)(2.10)

где k в нотации данной работы эквивалентно v_r . Первое слагаемое уравнения (2.10) отвечает в уравнении Аррениуса величине 2,301E_a/RT, и, соответственно, энергия активации будет равна 10,5 кДж/г-ат. О. Численное интегрирование¹ этой зависимости по интервалам температур позволяет получить величины v_ri, и затем распределить величину П/В по ступеням реактора по простому соотношению: 

.(2.11)

    Вид получаемой по уравнению (2.11) зависимости показан на рис.2.8. Более строгий расчет этих соотношений на основе гидродинамической модели конвективной системы [Тутубалин, Гричук, 1997а], в которой не использованы указанные выше допущения (а) и (б), дает весьма близкие результаты.

    Зависимость П/В от времени. Разрастание зоны метасоматизированных пород по нормали к поверхностям трещин, если оно контролируется скоростью поверхностной реакции, как это предполагалось выше, будет происходить по линейному закону от времени. Однако очевидно, что по мере расширения слоя метасоматизированных пород все большую роль должен играть диффузионный перенос вещества в этом слое, который будет замедлять процесс взаимодействия со свежей породой. Из теории диффузионного метасоматоза Д.С.Коржинского [Коржинский, 1969] следует параболический закон разрастания - ширина метасоматических зон меняется пропорционально  , где - время. Переход от линейного режима разрастания к параболическому в общем случае предсказать вряд ли возможно. Для моделей больших по объему и дебиту гидротермальных систем с П/В близким к 1 можно полагать, что линейная стадия заканчивается при прохождении уже первых порций раствора. Поэтому скорость разрастания слоя метасоматизированных пород на стенках трещин будет убывать:

,(2.12)

где  - скорость в начальный момент процесса, x - толщина измененного слоя. При пошаговой разбивке процесса взаимодействия, примененной в методе МПСР приращение П/В на каждой ступени реактора для j-ой волны можно записать в виде простого соотношения:

.(2.13)

    Из сопоставления уравнений (2.8) и (2.12) видно, что соотношение (2.13) предполагает постоянную величину площади поверхности взаимодействия, равную S_fr . Это справедливо для системы плоскопараллельных трещин, для других же видов геометрии блоков породы (шаровой отдельности, гексагональной системы трещин тепловой контракции) уменьшение v_r со временем будет более резким, чем задаваемое уравнениями (2.12) и (2.13).

    Предельный случай - это полное исчерпание свежей породы на ступени реактора, тогда величина П/В становится равной 0. Уравнения (2.12) и (2.13) этого не предусматривают, и в алгоритм моделирования следует помещать специальное ограничение на суммарный объем породы, могущей вступать в реакцию с породой на данной ступени реактора, вытекающее из соотношения:

,(2.14)

где n_cr - номер порции раствора, при которой достигается полное исчерпание свежей породы, Ф - трещинная пористость, _y и _w - плотности породы и воды соответственно. При измеренных для океанских гидротермальных систем величинах П/В [Von Damm 1990] 0,5B2 и трещинной пористости базальтов порядка 3^o/o [Hyndman, Drury, 1976] исчерпание свежей породы на высокотемпературных ступенях реактора достигается в модели при числе волн n 1000, то есть в долгоживущих системах. Это согласуется с наблюдениями на разрезах офиолитов [Richardson et al., 1987; Nehlig et al., 1994], где сплошная переработка пород в эпидозиты наблюдается на границе дайковых комплексов с изотропными габбро, маркирующей положение кровли магматической камеры (то есть самой горячей и химически активной части конвективной системы), тогда как вышележащие породы изменены довольно слабо.

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

    1. В модели предполагается однородность слоя метасоматизированной породы, тогда как согласно теории диффузионного метасоматоза этот слой представляет собой диффузионную метасоматическую колонку. Неучет этого явления ведет к тому, что в модели будет получаться равная миграционная подвижность для всех компонентов, не образовавших при данных П/В твердых фаз, в то время как природная картина будет более сложной - миграционная подвижность элементов по нормали к трещинным каналам будет дифференцированной. Поскольку в методе определения П/В по составу природных растворов [Von Damm et al., 1985] используются наиболее подвижные щелочные элементы. То для модели океанских гидротерм расчеты могут давать завышенное извлечение из пород некоторых других элементов (S, Pb, Zn) - в модели неявно задается их извлечение вплоть до фронтальной зоны диффузионной колонки, что может не соответствовать природной ситуации.

    2. Модель предполагает, что свежая порода реагирует с раствором как однородное вещество. Это соответствует использованию уравнения (2.10), полученному для состояния взаимодействия раствор-алюмосиликаты, далекого от равновесия. Однако это приводит к тому, что минеральный состав исходной породы, а главное - состав раствора не влияют на скорость этой реакции ². 

    3. Использованная в уравнении (2.13) величина (П/В)₁ относится к процессу, контролируемому диффузионным переносом, и распределение этих величин по температуре может не совпадать с задаваемым уравнением (2.11), полученным для процесса, контролируемого поверхностной реакцией. Температурная зависимость в уравнении Вуда-Уолтера отвечает энергии активации E_a гетерогенной реакции 10,5 кДж/г-ат О. Эта величина сопоставима с энергиями активации диффузии ионов в водном растворе, имеющими величину 16-19 кДж/ат. [Эрдеи-Груз, 1976], и энергиями активации диффузионно-контролируемых реакций в водной среде (< 20 кДж/моль [Lasaga, 1981]).

    Замена в уравнении (2.11) константы в подинтегральном выражении на характерную для диффузионного процесса не приведет к существенному изменению зависимости для П/В (см.рис.2.8). Отклонения в распределении П/В при этом будут меньше других погрешностей модели, и учет в уравнении (2.11) этих особенностей, усложняя модель, не даст улучшения ее общего качества.

    4. Принятый в модели способ разбиения трубки тока раствора на температурные участки предполагает полное смешение в пределах этого участка за счет гидродинамической дисперсии и отсутствие более масштабных процессов смешения растворов с разной "тепловой предысторией" внутри гидротермальной системы. В природных конвективных системах эти условия не всегда соблюдаются. Если выполнение первого из них может быть оптимизировано выбором температурного шага модели (по сути - размером блока породы, выделенного, как ступень реактора), то второе требует введения гидродинамической модели процесса и сеточного разбиения термодинамической модели в пространстве [Тутубалин, Гричук, 1997б].

5. Политермическая модель, основанная на методе МПСР, как и в рассмотренном выше случае изотермической колонки (рис.2.6), склонна к проявлению артефактов - колебаний состава раствора при изменении состава минеральных ассоциаций на предшествующих ступенях реактора.

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