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

3.3. Моделирование современного рудообразования в океане - результаты предшествующих исследований

    Впервые термодинамические расчеты, моделирующие высокотемпературное взаимодействие морской воды с базальтами выполнил Т.Волери в 1978 г., однако результаты этих исследований не были опубликованы.

    Термодинамическая модель поведения морской воды при нагреве до температур 300-350^oС, характерных для гидротерм океана была рассчитана в работах [Гричук, Борисов, 1980; Reed, 1982]. Результаты моделирования сопоставлялись с экспериментальными данными из работы [Bischoff, Seyfried, 1978], в которой было исследовано поведение морской воды в диапазоне 25-350^oС. Целью этих расчетов было доказательство возможности корректного описания в термодинамических моделях поведения морской воды в гидротермальном процессе. Хорошее согласие с экспериментом позволило перейти к более сложным задачам моделирования с участием силикатных пород.

    Расчеты изотермического взаимодействия морская вода-риолит были проведены М.Ридом [Reed, 1982], морская вода-базальт и морская вода-перидотит - автором данной работы совместно с М.В.Борисовым и Г.Л.Мельниковой [Гричук и др., 1982] и М.Ридом [Reed, 1983]. Эти расчеты показали удовлетворительное согласие с имеющейся экспериментальной информацией по соответствующим системам. Общие тренды изменения состава раствора также соответствовали наблюдаемым в природе. Вместе с тем, поскольку в указанных работах использовался метод "степени протекания реакции", в них не удалось получить правильной картины поведения серы (причина этого обсуждается в разделе 2.1), а вследствие этого - и рудных элементов.

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

    Моделирование рудоотложения из растворов "курильщиков" впервые было выполнено Д.Жанецки и У.Сейфридом [Janecky, Seyfried, 1984]. В этой работе был проведен расчет смешения гидротермальных растворов 21^oВТП с холодной морской водой. В модель было включено 15 элементов, в том числе Fe, S, Cu, Zn, Pb и Ba (состав мультисистемы по зависимым компонентам - минералам и частицам в публикации не указан). Из расчетов были исключены алюмосиликаты и силикаты: первые - из-за того, что в мультисистеме не задавался Al, вторые - из-за того, что в предварительных расчетах авторы получили слишком большие по их мнению количества осадков кварца, талька, серпентина и тремолита. Расчеты проводились в интервале температур от 350 до 25^oС при давлении 250 бар.

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

    Анализ примененной в работе методики моделирования показывает, что она представляет собой вариант метода "степени протекания реакции" Г.Хелгесона (см. главу 2), в котором происходит своеобразное "титрование" гидротермального раствора холодной морской водой без удаления продуктов реакции. Методика Д.Жанецки и У.Сейфрида принципиально не позволяет воспроизвести дифференциацию рудных компонентов. И действительно, в расчетах равновесного пути [Janecky, Seyfried, 1984, fig.12] в конце процесса, при температурах ниже 100^oС медные минералы (ковеллин и халькозин) присутствовали в том же количестве, что и халькопирит в начале процесса, тогда как при моделировании природного процесса медь в виде халькопирита должна была извлекаться из раствора и уже не входить в последующие минеральные ассоциации. Тем не менее, результаты моделирования в работе [Janecky, Seyfried, 1984] оказались во многом сходны с природными прототипами. Они показали возможность применения методов термодинамического моделирования для анализа даже таких относительно быстрых процессов, как гидротермально-осадочное рудообразование.

    Близкая по методу моделирования работа была выполнена Т.С.Бауэрс и др. [Bowers et al., 1985] на основании более полных данных по составу гидротермальных растворов 21^oс.ш. ВТП и бассейна Гуаймас. Мультисистема была дополнена рядом элементов - Al, Mn и Ag, а исходные анализы откорректированы на сдвиг рН при охлаждении анализируемых проб.

    Значительным шагом вперед в моделировании гидротермального рудообразования в океане явились работы М.Тиви и Р.Макдаффа [Tivey, McDuff, 1990; Tivey, 1995], в которых были рассмотрены процессы тепло-массопереноса и химического взаимодействия в стенках жерл курильщиков. В модели Тиви-Макдаффа учтены теплопередача, диффузия и течение раствора через пористую стенку и локальные термодинамические равновесия между раствором и твердыми фазами. Рассматривалась стационарная по температурному и концентрационным профилям задача.

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

Термодинамические модели формирования гидротермальных растворов и полные модели океанских гидротермальных систем.

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

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

    К числу общих проблем относятся создание согласованной базы термодинамических свойств веществ. В работах [Гричук и др., 1985; Гричук, 1988, 1996, Гидротермальные сульфидные руды , 1982] использовался банк данных UNITHERM (МГУ). Термодинамические свойства твердых веществ в нем взяты в основном по сводкам [Robie et al., 1978; Термодинамические свойства , 1978-1982] (см. также [Борисов, Шваров, 1992]). Для простых ионов в этом банке до 1989 г. использовалась работа [Helgeson et al., 1981], а после - модифицированная модель Хелгесона-Киркхама-Флауерса (HKF) [Shock, Helgeson, 1988, Johnson et al., 1992]. Свойства комплексных ионов вычислялись с помощью уравнения Рыженко в полуэмпирическом варианте [Методы геохимического моделирования , 1988; Борисов, Шваров, 1992] (см. Приложение). Для моделей М.Рида и Т.С.Бауэрс [Reed, 1983; Bowers, Taylor, 1985; Bowers, 1989] термодинамические свойства веществ были взяты в основном из [Helgeson et al., 1978, 1981; Robie et al., 1978]. В работе [Силантьев и др., 1992] термодинамические свойства веществ взяты из банка "ДИАНИК" (ГЕОХИ РАН), в котором для простых ионов используется модифицированная модель HKF, а для комплексных частиц - уравнение Рыженко в модификации О.В.Брызгалина [Брызгалин, 1985]. Из этого перечисления видно, что хотя различия между термодинамическими базами моделей разных авторов и существуют, они не очень значительны, поскольку базы опираются в основном на одни и те же фундаментальные сводки. Расхождения в результатах моделирования больше связаны с разными использованными наборами минералов и растворенных частиц. Существенные изменения в результатах моделирования произошли только после появления модифицированной модели HKF, значительно улучшившей точность расчета термодинамических свойств в отношении растворенных веществ при высоких температурах и повышенных давлениях.

    В отношении выбора составов исходных растворов и пород модели океанских гидротерм находятся в исключительно благоприятном положении по сравнению с моделями других типов гидротермальных объектов из-за большой однородности вмещающих пород и исходной морской воды. Весьма сходен и интервал условий - температура от 100-150 до 350-400^oС, давление до 500 бар. По этим параметрам указанные выше работы практически не отличаются (кроме [Силантьев и др., 1992] - см. ниже).

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

    Полная модель гидротермальной системы была построена в работе М.Рида [Reed, 1983]. В ней процесс рециклинга был представлен последовательностью из трех этапов: 1) нагревание морской воды до 300^oС при 500 бар; 2) изотермическая реакция воды с базальтом при 300^oС, 500 бар; 3) охлаждение раствора, сопровождающееся рудообразованием. На этапе нагрева морской воды из нее выпадает ангидрит, при этом удаляется 1/3 исходно присутствовавшего SO₄. Взаимодействие воды с базальтом рассчитывалось в работе М.Рида методом "степени протекания реакции". Полученные в модели тренды изменения концентраций Mg, Ca, SO₄, H₂S, Fe, Ag, Cu и Zn удовлетворительно согласовались с имевшимися в то время экспериментальными и натурными данными. Для процесса рудоотложения в работе М.Рида рассматривалось два сценария: а) кондуктивное охлаждение; б) смешение с морской водой при 25^oС. Для моделирования были взяты модельные растворы предыдущего этапа, полученные при П/В = 0,045, в котором H₂S > (Fe + Zn + Cu), и при П/В = 0,062, в котором соотношение было обратным. Охлаждение первого раствора дает осадок кварц + пирит + халькопирит + халькозин + сфалерит + акантит, при этом в воде сохраняется избыток сероводорода. При охлаждении второго к ассоциации добавляет магнетит, и сероводород полностью связывается. В сценарии со смешением халькопирит не образуется, а магнетит заменяется гематитом.

    Детальный анализ этой работы показал, что использование в ней двухэтапной схемы для нисходящей ветви (нагрев воды, затем реакция с базальтом) и метода "степени протекания реакции" привело к значительному искажению результатов. В модели М.Рида из-за возможности "back-reaction" была получена слишком большая область устойчивости сульфатной серы, вплоть до самых больших рассмотренных П/В, что привело к связыванию из раствора кальция. Концентрации сероводорода оказались высокими за счет восстановления морского сульфата, что находится в противоречии с изотопными данными по сере. Полученные М.Ридом соотношения между сероводородом и рудными металлами, когда с ростом П/В появляется дефицит сероводорода, как показали позднейшие данные, является артефактом метода моделирования. Это искажает все поведение халькофильных элементов. В сценарии рудоотложения при смешении не был учтен тепловой баланс - все полученные смеси приводились в равновесие при 25^o, хотя их температура должна зависеть от пропорции смешения - см. [Janecky, Seyfried, 1984].

    Несмотря на эти дефекты, работа М.Рида сыграла большую роль в формировании современных взглядов на химические процессы в океанских гидротермальных системах (см. [Кривцов, 1987]).

    Первая полная модель океанской гидротермальной системы, основанная на методе проточного ступенчатого реактора была опубликована в работах [Гричук, Борисов, 1982, 1983; Гричук и др., 1985]. В этой модели рассматривалось взаимодействие морской воды со свежим базальтом (то есть прохождение первой порции раствора). Для отношения порода/вода в модели была принята экспоненциальная зависимость от температуры, с таким расчетом, чтобы П/В составляла 0,173-1,143 кг/кг (см.рис.2.8) Для восходящей ветви и зоны рудоотложения было рассмотрено 3 сценария: а) медленное кондуктивное охлаждение, также представленное в виде проточного ступенчатого реактора; б) быстрое охлаждение - при излиянии на дно, но без химического взаимодействия с морской водой; в) охлаждение при смешении - аналог задачи, рассмотренной в работе [Janecky, Seyfried, 1984]. Из-за недостатка имевшейся к тому времени термодинамической информации модель была ограничена температурным интервалом 150-350^oС при Р_{нас.пара} и 13-компонентной мультисистемой, включавшей из рудных элементов - Fe, S и Pb. Мультисистема содержала 36 минералов и 48 растворенных частиц.

    Было получено общее соответствие минеральных ассоциаций метасоматитов и состава формирующихся в нисходящей ветви модельных растворов данным натурных наблюдений. При расчете быстрого охлаждения модельных растворов были получены минеральные ассоциации, типичные для взвешенного вещества "курильщиков". Однако в этом цикле работ не удалось достичь хорошего согласия с природными прототипами по некоторым компонентам раствора. В расчетах при 350^oС были получены заниженные концентрации Ca, слишком высокие концентрации Fe, Al и H₂, а величины pH_Т оказались слабо щелочными. Причины этих дефектов были связаны с недостаточно надежной термодинамической информацией по растворенным формам этих элементов (хлоридным комплексам Ca и Fe и гидроксидным комплексам Fe^III).

    Модель формирования океанских гидротермальных растворов была рассчитана в работе [Bowers, Taylor, 1985]. В ней нисходящая ветвь конвекции также была представлена как ступенчатый проточный реактор, и рассматривалось прохождение через него одной (первой) порции раствора. Термодинамическая мультисистема в модели включала 11 элементов, в том числе из рудных - только Fe и S, ее полный компонентный состав в работе не указан. Суммарное отношение порода/вода в модели было принято 1,787 кг/кг, что соответствует оценкам, полученным для гидротерм 21^oс.ш. ВТП (= 2 [Von Damm et al., 1985]). Для того, чтобы распределить эту сумму по ступеням модели, авторы использовали признаки правдоподобности результатов расчета по некоторым показателям минералогического и химического состава (появлению в ассоциациях минералов эпидота, возрастанию рН раствора), то есть решали по существу обратную задачу. Однако эти признаки в работе не формализованы, и выбор конкретных значений П/В для каждой ступени представляется в большой мере субъективным. Для расчетов равновесий в работе [Bowers, Taylor, 1985] использовалась программа EQ3/6. Поскольку эта программа использует предложенный Г.Хелгесоном алгоритм, главным параметром которого является - степень протекания реакции (см. главу 2), она не может рассчитать равновесие непосредственно для заданного состава системы порода-вода, а рассчитывает путь "титрования" раствора породой. Вследствие этого в работе Т.С.Бауэрс и Х.Тейлора для каждой температурной ступени была рассчитана закрытая система с величиной , меняющейся от 0 до заданного значения, после чего результирующий раствор переносится на другую ступень. Такой метод авторы назвали "complex open-system model", и считают его более реалистичным, чем другие методы моделирования¹. Путь такого процесса в координатах П/В vs Т представляет собой ломаную линию; он показан на рис.2.8. Процессы рудоотложения в восходящей ветви системы в работе [Bowers, Taylor, 1985] не рассматривались. Оригинальной особенностью модели Т.С.Бауэрс и Х.Тейлора является включение в нее изотопных равновесий (по H и O). Изотопная часть модели основана на методе Омото и использует результаты расчета химических равновесий в качестве входной информации (подробнее см. главу 5). Для реализации изотопных расчетов программа EQ3/6 была дополнена специальным блоком.

    В работе [Bowers, Taylor, 1985] было получено весьма реалистичное воспроизведение минеральных ассоциаций метабазитов океанской коры, и в целом удовлетворительное описание изменения состава морской воды при рециклинге. Однако, так же, как и в работе [Гричук и др., 1985], расчетные концентрации Fe и Al в растворе оказались завышенными, а Ca - сильно заниженными по сравнению с природными прототипами, при этом величины рН_Т были нейтральными вместо слабокислых. Принципиально важный результат получен в этой работе в изотопной части модели. Вычисленные для гидротермальных растворов модельные величины D = +2,67 ^o/oo и ¹⁸О = +2,0 ^o/oo хорошо согласуются с наблюдаемыми в природных гидротермах (см.табл.3.6). Это свидетельствует о том, что изотопный сдвиг гидротермальных растворов океана относительно морской воды полностью объясним реакцией с породами и не требует привлечения гипотез об участии в гидротермах каких-либо иных источников воды (в том числе - магматогенных флюидов).

    К сожалению, надежность полученной модельной оценки сдвига по D сейчас может быть подвергнуто сомнению. В работе Т.С.Бауэрс и Х.Тейлора в качестве исходной информации использовались экстраполированные значения для коэффициентов распределения изотопов водорода для пар минерал-вода. При этом, как опорные для определения вида экстраполяционной зависимости, были взяты экспериментальные данные по паре каолинит-вода, для которой Х.Тейлор, опираясь на имевшиеся тогда данные, получил волнообразную зависимость с максимумом около 300^oС. Соответствующий вид кривых он прогнозировал и для других водосодержащих минералов (кроме эпидота). Недавно было обнаружено [Gild, Sheppard, 1996], что низкотемпературные данные по каолиниту, использованные Х.Тейлором неверны - в этих опытах не было достигнуто изотопное равновесие. Вследствие этого изотопные сдвиги по водороду, полученные в модели Бауэрс-Тейлора могут быть завышены примерно вдвое. Однако эта погрешность не может повлиять на принципиальный вывод о роли морской воды, изложенный выше.

    "Изотопное направление" моделей океанских гидротерм получило дальнейшее развитие в работах [Гричук, 1988] и [Bowers, 1989], в которых было рассмотрено поведение в нисходящей ветви конвекции изотопов серы - одного из главных рудных элементов. Эти результаты обсуждаются в главе 5. В химическом отношении модели в этих работах были расширены путем включения в них Cu и Zn с соответствующими твердыми и растворенными формами.

    Метод "многоволнового проточного ступенчатого реактора" был применен для моделирования эволюции океанской гидротермальной системы во времени в работах [Краснов и др., 1990; Гричук, Леин, 1991; Гидротермальные сульфидные руды , 1992]. В указанных работах удалось улучшить согласие с природными прототипами по содержаниям в растворах Ca, Al, H₂ и величинам рН за счет использования при расчете исходных термодинамических данных модифицированной модели HKF.

    Оригинальный подход к описанию процессов взаимодействия вода-порода был предложен в работе [Силантьев и др., 1992]. Здесь также использовалась идея многоволнового проточного ступенчатого реактора. Основная задача работы состояла в воспроизведении метаморфизма пород океанической коры, поэтому авторы рассмотрели неоднородный разрез коры, сложенный базальтами, габброидами и гипербазитами. При этом они отказались от принципа частичного равновесия и полагали, что весь объем пород эффективно участвует во взаимодействии. Соответственно строилась и процедура расчетов: каждая порция воды на каждой ступени реагировала с большой массой породы (10 кг породы/кг раствора, что примерно отвечает пористости базальтов), после чего переносилась на следующую ступень. Модель рассчитывалась в диапазоне 150-500^oС с шагом 50^o, при этом ступени до 350^o включительно содержали базальты, а далее - габбро, либо гипербазиты. Мультисистема включала 14 элементов (в том числе Fe, S и Pb), образующих 45 твердых фаз и 39 растворенных частиц. Для расчетов использовалась программа EQUILIBR. Расчет был закончен после прохождения 60 порций воды, в результате интегральное отношение порода/вода для первой ступени составило 0,167, а для последней "базальтовой" ступени (350^o) - 0,83.

    Из сравнения с материалами, приведенными в разделе 2.1 видно, что такой метод (ступенчатый реактор с большими массами породы на ступенях) должен удовлетворительно воспроизводить метасоматическую колонку инфильтрационного метасоматоза, однако для получения состава тыловых зон колонки требуется гораздо большее количество порций раствора (П/В должно стать по крайней мере менее 0,02). Поэтому в результатах расчета для ступени 150^oС, интерпретируемой в работе [Силантьев и др., 1992] как пример цеолитовой фации метаморфизма, преобладают тремолит и эпидот. По-видимому, для верхней части базальтоидного разреза представление об эффективном взаимодействии с водой всего объема пород не вполне согласуется с реальностью. В то же время для более горячих частей модели в работе получено хорошее соответствие модельных ассоциаций природным, типичным для метабазитов зеленосланцевой фации метаморфизма.

    Вместе с тем данные по составу модельных растворов, приведенные в [Силантьев и др., 1992], расходятся с наблюдаемыми в океанских гидротермах по многим показателям. Из приведенных результатов для ступеней 350 и 400^oС, наиболее близких к реализуемым в природных гидротермах условиям, только для Na наблюдается соответствие природе, а для K, Ca, SiO₂, Pb концентрации занижены на 0,5-2 порядка. Обращает на себя внимание и преобладание окисленных форм серы над восстановленными, чего в гидротермах океана не наблюдается. Это свидетельствует о существовании у авторов работы проблем с термодинамическим описанием поведения растворенных компонентов. Одновременно это показывает, что при верификации моделей использование данных только по минеральным ассоциациям недостаточно, модель должна проверяться также и по составу раствора.

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

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